Identification and ordering of drainage divides 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. We thus propose a divide-network metric that orders divides based on the average distance one would have to travel down on either side of a divide to reach a common stream location. Because measuring these distances is computationally very expensive, we instead sort divide segments in a tree-like network, starting from endpoints at river junctions. The sorted nature of the network allows assigning distances to points along the divides, which can be shown to scale 15 with the average distance downstream to the common 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 a natural landscape and to results from landscape evolution model experiments to assess which parameters of divides and divide networks are diagnostic of divide mobility. Results show that stable drainage divides strive to attain a constant hillslope relief as well as flow distance from the nearest stream, provided a distance of > ~5 km from 20 endpoints. Disruptions of such pattern can be related to mobile divides that are lower than stable divides, closer to streams, and often asymmetric in shape. In general, we observe that drainage divides 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 for a longer time period. We compared different topographic metrics to assess drainage divide mobility and found that cross-divide differences in hillslope relief proved more useful for assessing divide migration than other tested metrics. Finally, we 25 introduced a new metric to assess divide junction stability, based on the connectivity of junctions with adjacent drainage divide segments. https://doi.org/10.5194/esurf-2019-51 Preprint. Discussion started: 2 October 2019 c © Author(s) 2019. CC BY 4.0 License.


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).At regional scales, mobile divides can lead to important changes in 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).The lateral migration of divides is a consequence of spatial gradients in surface uplift (positive or negative) as well as stream captures, which frequently accompany tectonic deformation due to shearing, stretching, and rotating stream networks (Bonnet, 2009;Castelltort et al., 2012;Willett et al., 2014;Goren et al., 2015).
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).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 cross-divide differences of χ values, a proxy for steady-state river channel elevation (Perron and Royden, 2013).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 lowrelief 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 unadjusted parts of landscapes.They also cast doubt on the ease of comparing crossdivide 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.Amongst the topographic parameters that they tested are cross-divide differences in channel elevation at a reference drainage area, mean headwater hillslope gradient, and mean headwater local relief.Focusing on the stability of divides, Spotila (2012) argued that divide junctions and pyramidal peaks are more stable than linear solitary divides and might act as anchor points for drainage divide networks.
Whether based on the geometry of hillslopes or river channels, previous studies on mobile divides generally focused on relative differences in either one across several specific, manually selected drainage divides (e.g., Willett et al., 2014;Goren et al., 2015;Whipple et al., 2017;Buscher et al., 2017;Beeson et al., 2017;Gallen, 2018;Guerit et al., 2018;Forte and Whipple, 2018).Even if appropriate in these studies, such a procedure introduces unwanted subjectivity, both in the selection of divides and how any cross-divide comparison is done.This choice of procedure may be attributed to the fact that, so far, no straightforward approach exists to reliably extract the drainage divide network from a digital elevation model (DEM).
Functions that classify topographic ridges (the common shape of drainage divides) based on local surface characteristics and a threshold value (e.g., Little and Shi, 2001;Koka et al., 2011) are prone to misclassifications.The grey-weighted skeletonization method by Ranwez and Soille (2002) (homotopic thinning) requires the determination of topographic anchors (e.g., regional maxima), which makes it sensitive to DEM errors.The approach by Lindsay and Seibert (2013), who identified pixels belonging to drainage divides based on confluent flow paths from adjacent DEM pixels and a threshold value, is computationally expensive and sensitive to edge effects that depend on DEM size.Furthermore, drainage divides that coincide with pixel centres are inconsistent with the commonly used D8 flow routing algorithm (O'Callaghan and Mark, 1984), in which each pixel belongs to a specific drainage basin.A probabilistic approach, based on multiple flow directions exists (Schwanghart and Heckmann, 2012), but computation is expensive and thus restricted to a few drainage basin outlets.Finally, all of these approaches merely yield a classified grid, but no information about the tree-like network structure of drainage divides, which requires ordering of the divide pixels into a network (Figure 1).
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 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., Horton, 1945;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 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 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.In the following, we first present an approach that allows identification and ordering of drainage divides in a DEM.We investigate ways of ordering drainage divide networks and analyse basic statistical properties with a natural example from the San Gabriel Mountains in Southern California.We then present experiments with a numerical landscape evolution model that we conducted to investigate how drainage divide networks respond to external perturbations.
Finally, we return to the San Gabriel Mountains and analyse the drainage divide network for signs of drainage divide mobility.

Drainage divides in digital elevation models
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 (O'Callaghan and Mark, 1984) 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 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.As a result, drainage basin boundaries, and thus divides, are located in between DEM pixels, and have infinitesimal width.Another important 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.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 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).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 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 arise near confluences (Figure 2).
Just like streams, drainage divides are organized into tree-like networks (Figure 1).Whereas streams start at channel heads, divides start (and end) at stream confluences and follow the boundaries between adjacent drainage basins.We refer to these points as endpoints of the divide network (Figure 2).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.In stream networks, flow is directed and leads to natural sorting of stream nodes by elevation and the hierarchy of stream segments can be easily described by their upstream area, for example.In contrast to streams, however, there exists no obvious metric for ordering and classifying divides.We propose that a meaningful metric for ordering divides is the average branch length (Lindsay and Seibert, 2013), 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 4).Measuring this distance 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 (Lindsay and Seibert, 2013), it is not very practical.However, as we will show later, this distance scales linearly with the maximum distance along the divide network, as well as the maximum number of divide segments or junctions, both of which are more easily computed.We thus propose to order the nodes and edges of the divide network, by their maximum down-divide 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 divide network the divide distance (dd).

Divide algorithm
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.For a given DEM, we first defined a stream network, based on the D8 flow direction algorithm and a minimum upstream area.We then applied the following steps.
(1) We defined drainage divides based on drainage basin boundaries that we obtained for upstream areas at tributary junctions and drainage outlets.These divide segments do not cross any rivers but their nodes may coincide with stream edges (Figure 2).
(2) We removed redundant divide segments in the collection of divides, which arise from nested and adjoining drainage basins.As a result, we are left with a set of unique divide segments, which, however, may be continuous across junctions or terminate where they should be continuous (Figure 3).
(3) We organized the collection of divide segments in a divide network.This is the core of the algorithm, in which we effectively identify endpoints and junctions, merge broken divide segments, and split divide segments at junctions (Figure 3).Our algorithm distinguishes between junctions, endpoints, and broken divide segments, based on the number of edges, the number of segment termini, and the existence and direction of a diagonal flow direction at the terminal nodes of divide segments (Appendix 7.1).
(4) Finally, we sorted the drainage divide segments of the network.We iteratively identified segments that are connected to endpoints and removed them from the list of divide segments until no divide segments were left (Figure 5).
After sorting, we computed maximum distances along the divide network and assigned orders to divide segments, based on the ordering of stream networks, first introduced by Horton (1945).We adopted both Strahler's (1954) andShreve'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; ω 1 ) Eq. 1 Eq. 2 Topo:  " = max( + ) + 1 Eq. 3, where i are the divide orders of the  joining divide segments and k is the divide order of the following divide segment.In 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 3).Based on the Strahler ordering scheme, the bifurcation ratio  9 , can be derived from (Horton, 1945): where  ; is the number of divide segments of order .
Our divide algorithm currently does not handle internally drained basins very well.Whereas the divides of internally drained basins are easy to identify, they are not easily sorted in a meaningful manner.In their case our proposed metric of ordering drainage divides, based on the distance to a common stream location, is undefined.At the moment, such divides are simply left unsorted and are flagged as belonging to an internally drained basin.

Topographic data and analysis
We investigated basic characteristics of drainage divide networks using a 30-m resolution DEM from the 1-arc second Shuttle Radar Topography Mission data set (Farr et al., 2007).We focused on the catchment of the Big Tujunga River in the San Gabriel Mountains, USA.The catchment is a good example of a transient landscape with active drainage basin reorganization and landscape rejuvenation as the river incises into a relict pre-uplift landscape (DiBiase et al. 2014).We pre-processed the DEM by carving through local sinks (Soille et al. 2003) to avoid artificial internally drained basins.
We analysed the divide network both on its planform geometry and on the topography it is related to.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 are subsequently associated, or mapped, to divide edges and entire divide segments.As topographic metrics, we included elevation (z), as well as hillslope relief (HR) and horizontal flow distance to the stream network (FD).We also computed χ values on either side of a divide, using a reference area A0 of 1 m 2 , a reference concavity θref of 0.45, and setting the base level xb to 0 at the edge of the model domain (e.g., Perron and Royden, 2013): Eq. 5.
For each divide edge, we computed these topographic metrics for the two neighbouring pixels that belong to adjacent drainage basins (Figure 2) and denoted the cross-divide minimum, maximum, sum, difference, and average in any one topographic metric , as  M+4 ,  MNJ , ∑ , ∆, and  P , respectively.Topographic metrics of entire divide segments are based on those of the divide edges that it is composed of.For quantifying cross-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 related to the cross-divide difference in hillslope relief, and we define its absolute value as the divide asymmetry index (DAI): Spotila ( 2012) suggested that the stability of divide junctions is related to the number of joining drainage divides.Because divide junctions, as defined here, cannot connect more than four divides segments -and most often connect three segmentswe introduce a new metric to quantify divide junction connectivity, CJ: We define CJ to correspond to the ratios of the Euclidean distance, d, and the divide distance, dd, integrated over a specified maximum divide distance, dd,max, normalized by dd,max.The quantity CJ is sensitive to the number of divides within a given divide distance from a junction, weighted by their orientation towards the junction (Figure 6).The value dd,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, CJ is not sensitive to the value of dd,max.However, for actual junctions, CJ is typically sensitive to the value of dd,max, because as dd,max grows, increasingly more junctions are at a distance dd<dd,max of a specific junction and thus the number of divide segments grows with dd (Figure 6a).In general, CJ will be sensitive to the position of a junction within the drainage divide network, if the junction's maximum divide distance from an endpoint is smaller than dd,max.In other cases, CJ 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.

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): where 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 10 Myr, until the topography reached nearly steady state.The result of this model, which we termed 'Initialize', provided the starting point for four other models that we again ran for 10 Myr (Figure 7).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).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

Basic divide statistics
We applied our divide algorithm to the Big Tujunga catchment, based on a stream network that we derived from a minimum upstream area of 0.1 km 2 .The resulting divide network for different ordering schemes and the associated frequencies of exceedance of divide orders are shown in Figure 8.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 junction that is encountered last in the ordering process (black arrows in Figure 8).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 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, Shreve divide orders increase more rapidly along the northern perimeter compared to the southern, and the junction encountered last during the sorting process opposes divide segments of different orders (Figure 8).
For the Strahler ordering scheme, the frequency of exceedance of divide segments decreases exponentially with divide order (ω), which is consistent with Horton's law of stream numbers (Horton, 1945), and corresponds to a bifurcation ratio of  9 = 3.89 ± 0.30 (standard error).The bifurcation ratio of the associated stream network is 5.39 ± 0.87.The average length of divide segments is positively correlated with the minimum upstream area used for deriving the stream network, simply because both the stream and the divide network extend to finer branches.For the value that we used (0.1 km 2 ), the average length across all divide orders is 442 ± 323 m (±1σ).The average length for different divide orders tends to be lower at ω < ~40 (based on the Topo ordering scheme), compared to ω ≥ ~40 (Figure 9).We quantified the proposed divide metric, 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 (Λ), for 100 randomly chosen divide edges per divide order in the Topo ordering scheme.Although the maximum order, for which Λ can be determined in our DEM, is limited to ω ≤ 55, results demonstrate that Λ (in km) increases linearly with 0.36 × divide order (ω), and 1.11 × maximum divide distance (dd, in km) (Figure 10).The linear scaling of these two relationships is a consequence of the similarity of segment lengths for different divide orders (Figure 9).

Numerical experiments
The simulated landscapes along with their drainage divide networks at the end of the numerical experiments are shown in Figure 11 and Figure 12, and in the Supplementary Material we provide movies of all simulations.Figure 13 shows how the frequency of exceedance of divide lengths for different orders changed during the landscape evolution experiments.The first 1-2 Myr of the Initialize model were characterized by large changes in divide lengths and orders, but relatively few changes thereafter.Compared to the Reference model, in which the divide network structure did not change much after ~5 Myr, the Rotating, Inclined, and Spheres models exhibit changes in the divide network, mostly at higher divide orders.To provide a measure of the mobility of drainage divides, we computed the percentages of drainage area that were exchanged between catchments during the simulations (Figure 14).Except for the Reference model, all models are characterized by significant changes in drainage area and thus mobile drainage divides.Area changes in the Initialize model were large in the beginning but levelled off rapidly during the first 1 Myr.In the Rotating model, large area changes appeared as discrete pulses, induced by drainage captures of major streams (Figure 12b), whereas the background area changes during the rotation/faulting were relatively small (<0.1%/40kyr).Area changes in the Inclined model were moderate (~0.25%/40kyr) throughout the simulation and oscillated in conjunction with the passage of more resistant layers through the landscapes (Figure 12c).Area changes in the Spheres model were generally more pronounced if the resistant spheres appeared in the course of rivers, which forced them to steepen and induced surface uplift upstream, as compared to their appearance at drainage divides (Figure 12d).specific divide distance interval (Figure 15a-d).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.In the Spheres model, this pattern is sporadically disrupted (e.g., at ~13-15 Ma), when spheres emerged at the surface in places where they affected only a small portion of the drainage divide network.
Amongst the studied parameters, we find that mean elevation (z), hillslope relief (HR), and flow distance (FD) approach constant values in the Reference model, but are prevented to do so in the disturbed models.Temporal variations in junction connectivity (CJ), which captures the relative connectivity of drainage divide junctions for a dd,max value of 5 km, is similarly high in the Reference model as in the disturbed models.Together with the mean elevation it is the only parameter, where the actual magnitude of CJ depends strongly on the divide distance, simply because of the chosen dd,max value (see section 3.2).
Whereas the actual value of FD appears to show some dependency on divide distance, no dependency can be seen for HR.
Comparison of cross-divide differences in erosion rate with differences in hillslope relief, flow distance and chi (Figure 15eh) shows that absolute and normalized differences in hillslope relief (i.e., the divide asymmetry index), appear to be the most sensitive to drainage divide mobility, whereas cross-divide differences in chi and flow distance are sensitive to divide mobility in the Rotating model, but less so in the Inclined and Spheres models.
Figure 16 shows how minimum hillslope relief (HRmin), minimum flow distance (FDmin), and the divide asymmetry index (DAI) vary with divide distance in the Reference model and the three models with induced landscape disturbances.In contrast to the average values in Figure 15, we provide these metrics for all divide edges during the last 0.8 Myr of the model runs.
Note also that we sorted the data points so that high DAI values are in the front, to better assess where asymmetric divides are located, but that in all four models, high DAI points plot on top of low DAI points.In the Reference model, both minimum hillslope relief and minimum flow distance reach relatively steady values (HRmin ~250-350m; FDmin ~400-600m) at a divide distance of ~10 km.At low divide distances, both HRmin and FDmin approach zero -because divides are defined to start at the stream network -and these divides become increasingly more asymmetric.It is notable however, that some of the highest HRmin and FDmin values are also observed at low divide distances of approximately 1-2 km.The cluster of highly asymmetric divides with low HRmin and FDmin values at ~5 km divide distance is related to isolated adjustments in the lower right corner of the model (cf.videos in the Supplementary Material).In the three other models, the transition between quite variable divides at short distances and more steady 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 HRmin and FDmin values at higher distances.
These divides are particularly prominent at a distance of ~15-20 km and correspond to the position of the strike-slip fault.
Where  but are all shifted to higher elevations (cf., Figure 11).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 CJ values is equal to that of the Reference model, but high CJ values are steered towards the less erodible spheres (Figure 17a).

Anomalous drainage divides in the Big Tujunga catchment
With the information of how the different divide metrics recorded landscape disturbances in the numerical experiments, we returned to the Big Tujunga catchment to assess whether the divide network shows signs of disturbances.Notable deviations occur at ~22-25 km, ~ 41-45 km, and >54 km divide distance, and are typically associated with asymmetric divides (Figure 18). Figure 19a shows the geographic position of anomalously low (HRmin <200 m) and asymmetric divides (DAI>0.5),which are less than 1000 m from a stream (FDmin <1000m), and Figure 19b-e give impressions of how they appear in Google Earth.Most of these divides can be seen to border regions of contrasting local relief (Figure 19a), and many cluster along the eastern edge of the catchment.The predicted migration direction of divides along the southern and eastern edge of the catchment, from higher to lower relief, is mostly inwards, causing area loss.Furthermore, the junction connectivity along the eastern edge is relatively low, despite being high up in the divide network (Figure 20a).While the overall trend is that divides at greater divide distance (dd) are higher in elevation and have higher junction connectivity (Figure 20b), in the Big Tujunga catchment, the highest junction connectivity is associated to a divide segment in the southwestern part of the catchment (Figure 20a), which is also high in elevation and local relief (Figure 19a).

Extraction and ordering of drainage divide networks
Our new approach allows routinely extracting drainage divide networks from any DEM.We have shown that the maximum divide distance dd, calculated as the maximum distance along the 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 10).In contrast to the average branch length, however, the divide distance is more easily and rapidly calculated and is less prone to edge effects (cf., Lindsay and Seibert, 2013).The proposed sorting procedure (Figure 5) 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 lengths across all divide orders (Figure 9), divide orders derived with the Topo ordering scheme can serve a similar purpose as divide distance.
Alternatively, divide orders derived with the Strahler ordering scheme, can be used to investigate how the divide network conforms to Horton's (1945) laws of network composition.In the Big Tujunga catchment, for example, the bifurcation ratio of the divide network (Rb ~3.9) is lower than that of the stream network (Rb ~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 Rb ~4.6 for the divide network.We obtained similar bifurcation ratios for the divide (Rb ~3.74 ± 0.48) and stream (Rb ~4.53 ± 0.18) networks in the Reference model of the numerical simulations, when taking only the four largest, similar-sized stream networks into account.The above 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 geometry.However, more observations from different landscapes are needed to assess systematic differences in between divide and stream networks.

Divide network dynamics
Stream networks tend to attain configurations that are in equilibrium with the geological and climatic environment, given an initial condition (e.g., Rinaldo et al., 2014).Because drainage divides result from the adjacency of 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 13), and divide distances.During the subsequent extension and shrinkage of 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 13).
In general, divide segments of high order, i.e., at great distance from endpoints, appear to be the most responsive to landscape disturbances (Figure 15).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 11b).In the cases of the Inclined and Spheres models, it may be 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 divide 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 16).In the other models, most of the investigated morphometric parameters are quite variable over the same distance and can be seen to rapidly adjust to disturbances such as drainage captures or migrating divides.Such behavior 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, appears appropriate.
This minimum divide distance could be lower or higher, depending on factors like drainage density and average hillslope relief, for example.

Quantifying drainage divide mobility
The analysis of stream networks has become a standard tool for inferring tectonic forcing and landscape history (e.g., Wobus et al., 2006;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 do they tell us about landscape history?Our Rotating model induced relatively sudden drainage captures (Figure 14).Because such events are always 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 16).More gradual divide migration, however, likely lacks such simple diagnostic criteria, and in those cases, cross-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 cross-divide difference in χ (Willett et al., 2014).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 that we included in this study.
Figure 21 shows how normalized cross-divide differences in χ, hillslope relief (HR), as well as flow distance (FD) compare to normalized cross-divide differences in erosion rate (ER), evaluated for each divide edge from the last 2 Myr of the landscape evolution experiments Reference, Rotating, Inclined, and Spheres.We find that cross-divide differences in hillslope relief (HR; Figure 21a) are most sensitive to the disturbances included in the models, whereas cross-divide differences in χ are similarly sensitive to disturbances in the Rotating model, but less so in the Inclined and Spheres models (Figure 21b; Figure 15).This difference may be explained by the fact that in the latter two models, we introduced spatial variability in the erosional efficiency of rivers (Kr) that we did not account for in our cross-divide comparison of χ, as would be required (Willett et al., 2014).Cross- 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 21c).However, there also exists substantial scatter in the relationship between crossdivide 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 21), which corresponds to the value below which we observe more variability in divide morphology, even in the Reference model (cf., Figure 16).We speculate that this 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 and Reference models 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 12 and studying the videos of the landscape evolution experiments (see Supplementary Material) is insightful: where drainage divides are migrating, one typically observes a range of cross-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 cross-divide differences in topographic metrics lagging behind cross-divide differences in erosion rate.
Our new junction connectivity index (CJ), complements existing topographic metrics in assessing drainage divide stability.
For example, the junction connectivity in our Rotating model is low along the fault (Figure 17), 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 12c, Supplementary Material).In this case, parts of the drainage divide network were mobile, not stable, but they moved towards particularly stable portions in the landscape.Therefore, the junction connectivity index (CJ) may also be interpreted as an attractor index that quantifies how strong a drainage divide network is pulled towards and anchored at a certain junction (Spotila, 2012).
As a final note, we emphasize that the above conclusions focus on the results from our numerical experiments, which depict an idealized world.It is clear that complexities present in nature, such as anisotropic and variable rock properties, hydroclimatic gradients, or biological influences on erosion processes and rates can lead to landscape patterns that bias any of the above topographic metrics, and need to be taken into account when interpreting both stream and divide networks.

Drainage divide mobility in the Big Tujunga catchment
Based on the observation of characteristic values of hillslope relief (300-500 m) and flow distance (1000-1800 m), we identified drainage divides in the Big Tujunga catchment that are anomalously low, close to a stream, and asymmetric (Figure 18, Figure 19), suggesting that they are mobile.Anomalous drainage divides are particularly frequent along the eastern edge of the catchment, where an area of low hillslope angles and local relief (Figure 19), 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 Mountains, and is currently destroyed by headward incision of rivers (Spotila et al., 2002;DiBiase et al., 2014).Cosmogenic 10 Be-derived erosion rates confirm lower erosion rates in the Chilao (Figure 20).
We identified another stretch of anomalous divides along the southern margin of the Big Tujunga catchment (Figure 19a,d), part of which is coincident with the trace of the San Gabriel Fault, which follows the orientation of the valley (Morton and Miller, 2006).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 (Scherler et al., 2016), suggesting weaker rocks closer to the fault (e.g., Roy et al., 2015).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 19b,d).We thus suggest that most, if not all of the anomalous divides we identified based on hillslope relief, flow distance, and divide asymmetry, are in fact 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 profiles in a tectonic framework (DiBiase et al., 2014).

Conclusions
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 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.
Based on landscape evolution model experiments in which we forced divides to migrate, we found that stable drainage divides 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 and in the case study of the Big Tujunga catchment, San Gabriel Mountains, 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 for a longer time period.In our comparison of different topographic metrics to assess drainage divide mobility, we found that cross-divide differences in hillslope relief proved more useful for assessing divide migration than other tested metrics.Finally, we introduced a new metric to assess divide junction stability, based on the connectivity of junctions with adjacent drainage divide segments.

Classification of divide network nodes
Once the drainage divides are defined, based on the outline of drainage basins, and redundant divide segments are removed, 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 xi 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 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.
is a parameter of the efficiency of river incision (T -1 ) (Kr = 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 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 ofWhipple et al. (2017b), who tested a wide range of rock uplift and erosional efficiency parameters and https://doi.org/10.5194/esurf-2019-51Preprint.Discussion started: 2 October 2019 c Author(s) 2019.CC BY 4.0 License.
<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 (Kr), dipping 30° towards northwest.The 'Inclined' model is representative for a landscape, in which rivers incise into tilted sedimentary rocks of non-uniform rock strength, similar to what has been studied byForte and Whipple (2018).During the experiment, the combination of surface lowering and inclination resulted in the strong layers to regularly sweep from southeast to northwest across the simulated landscape.The model 'Spheres' included 30 randomly assembled spheres of 3 km diameter with 75-% reduced erosional efficiency of rivers (Kr).This experiment may represent incision of a region that is characterized by country rocks with more-resistant magmatic intrusions.The expected behaviour of this model is similar to the landscape response to localized perturbations studied byO'Hara et al. (2018).

Figure 15
Figure15shows how the above-described disturbances affected drainage-divide metrics of simulated landscapes during the simulations.For all models, we computed the averages of several topographic parameters, measured at drainage divides of a HRmin and FDmin are low, DAI values are relatively high, although there also exist divides that have high DAI but regular HRmin and FDmin values.In the Inclined and Spheres models, the HRmin and FDmin values are never as low as in the Rotating model at divide distances >15 km, which reflects the lack of drainage captures.Instead, we observe frequent excursions to both higher and lower HRmin and FDmin values, either across all divide distances (Inclined) or at specific locations (Spheres).Deviations from the average values in the Reference model are greater in the Spheres model compared to all others.https://doi.org/10.5194/esurf-2019-51Preprint.Discussion started: 2 October 2019 c Author(s) 2019.CC BY 4.0 License.In the three disturbance models, DAI values are generally higher than in the Reference model (except for the cluster at 5 km divide distance) -although divides with high HRmin and FDmin values and at great divide distances always have somewhat lower DAI values.

Figure 17
Figure17shows the spatial pattern of divide junction CJ values at the end of the landscape evolution experiments.In general, junctions with higher CJ values tend to occur at higher elevation and at greater divide distance (dd).In the Reference model, the highest CJ values occupy the centre as well as the centres of the four quadrants of the model domain, mimicking the five of a six-sided dice.In the Rotating model, these clusters of high CJ values are maintained, but their connection is disrupted by the strike-slip fault, which induces low CJ values.The centrally located divide junctions occupy a similar range in CJ values, Figure 18 is similar to Figure 16 and shows how minimum hillslope relief (HRmin), minimum flow distance (FDmin), and the divide asymmetry index (DAI) vary with distance along the divide network of the Big Tujunga catchment.As in the simulations, we observe the largest range of HRmin and FDmin values, as well as some of the most asymmetric divides at low divide distances (<5 km), whereas most of the divides at higher distances cluster around HRmin values of 200-500 m and FDmin values of 800-1800 m.

Figure 1 :
Figure 1: Big Tujunga catchment, San Gabriel Mountains, United States.(a) Stream network, and (b) drainage divide network, both draped over hillshade image.Drainage divide network obtained with the approach developed in this study.Thickness of the stream and divide lines is related to upstream area and divide order, respectively.

Figure 2 :
Figure 2: Definition of drainage divides in a digital elevation model.Note the point where a drainage divide (red) coincides with a river channel (blue).See text for details.

Figure 3 :
Figure 3: Transformation of a collection of drainage divide segments (a) into a drainage divide network with endpoints and junctions (b).Black lines are drainage divides, blue lines are streams, and flow directions are shown as light gray lines.

Figure 4 :Figure 5 :
Figure 4: Ordering of divides based on the average distance to a common stream location, Λ. Numbered circles are places on drainage divides, and blue lines indicate the flow path to their common stream location (red circle).The resulting order of the drainage divide places is 1 < 2 < 3.

Figure 6 :
Figure 6: Junction connectivity (CJ) of different drainage divide junctions.(a) Map-view of four hypothetical divide junctions.Grey circles indicate Euclidean distance from the junction.Numbers to the southeast in each map gives the CJ value, based on a dd,max of 3000 m.(b) CJ values by divide distance for each case in (a).

Figure 7 :
Figure 7: Graphical representation of the model setups in the landscape evolution experiments.(a) 'Initialize' started from a nearly flat surface and reached almost steady state topography; the end of 'Initialize' provided the starting point for all other models.(b) 'Reference' included no changes.(c) 'Rotating' included a circular left-lateral strike slip fault, active throughout the experiment.(d) 'Inclined' included 1-km thick and 5-km spaced layers of 50-% reduced erodibility, dipping 30° towards northwest.(e) 'Spheres' included 30 randomly assembled spheres of 3 km diameter with 75-% reduced erodibility.

Figure 8 :
Figure 8: Divide network of the Big Tujunga catchment in the western San Gabriel Mountains, California, USA.Left column shows the divide network, with line thickness indicating the divide order.Arrow marks the last divide segment encountered in the sorting process.655

Figure 9 :
Figure 9: Average length (±1σ) of drainage divide segments by order (Topo ordering scheme).The number of observations per divide order drops below 10 at an order of 25. 660

Figure 10 :
Figure 10: Distance to common stream location by (a) divide order, and (b) divide distance, for 100 randomly chosen divide edges per divide order within the Big Tujunga catchment.

Figure 11 :
Figure 11: Modelled topography and drainage divide network (red solid lines) at the end of the landscape evolution experiments: (a) Reference, (b) Rotating, (c) Inclined, and (d) Spheres.White stippled lines show strike-slip fault in (b) and regions with reduced erodibility in (c) and (d).

Figure 12 :
Figure 12: Erosion rates and divide asymmetry index (DAI) at the end of the landscape evolution experiments: (a) Reference, (b) Rotating, (c) Inclined, and (d) Spheres.Black stippled lines show strike-slip fault in (b) and regions with reduced erodibility in (c) and (d).

Figure 13 :
Figure 13: Changes in cumulative divide length for divides of different orders (ω) during the landscape evolution experiments.Divides have been ordered with the Topo ordering scheme.

Figure 14 :
Figure 14: Changes in drainage area during the landscape evolution experiments.Y-axis unit is in percent per time step, which was 10 kyr.680

Figure 16 :
Figure 16: Minimum hillslope relief (a) and minimum flow distance (b) of all divide segments along the drainage divide network from the 690

Figure 17 :Figure 18 :
Figure 17: Divide junction connectivity (CJ) from the final stage of the landscape evolution experiments Reference, Rotating, Inclined, and Spheres.CJ was calculated with a dd,max value of 5000 m.(a) Map-view of divide junction CJ values.(b) Scatter plot of CJ by elevation, coloured by natural logarithm of the divide distance (dd).Black dashed line indicates area occupied by divide juntions in the reference model.

Figure 19 :
Figure 19: Anomalous divides in the Big Tujunga catchment.(a) 1000-m radius local relief map draped over hillshade image.White lines show the divide network and red lines depict asymmetric (DAI>0.4)divide edges with minimum hillslope relief <200m and minimum flow distance <1000m.Black arrows indicate the direction and magnitude the DAI, with the arrow pointing in the direction of lower relief, i.e., the inferred direction of divide migration.(b-e) Oblique views of asymmetric divides shown in (a), derived from © Google Earth.

Figure 20 :
Figure 20: Divide junction connectivity (CJ; dd,max = 5 km) of the Big Tujunga catchment.(a) Map-view of CJ values (colored circles), draped over divide network (gray lines), with highlighted anomalous divide segments (red lines).(b) Scatter plot showing the connectivity, elevation and divide distances (dd) of the junctions shown in (a). 715

Figure 21 :
Figure 21: Relationship between cross-divide differences in the topographic metrics hillslope relief (HR), chi (χ), and flow distance (FD), with cross-divide differences in erosion rates (ER) for all divide edges in the last 2 Myr of the landscape evolution experiments Reference, Rotating, Inclined, and Spheres.Marker colour denotes the divide distance (dd).