## Structured Abstract

### Introduction

River networks, the backbone of most landscapes on Earth, collect and transport water, sediment, organic matter, and nutrients from upland mountain regions to the oceans. Dynamic aspects of these networks include channels that shift laterally or expand upstream, ridges that migrate across Earth’s surface, and river capture events whereby flow from one branch of the network is rerouted in a new direction. These processes result in a constantly changing map of the network with implications for mass transport and the geographic connectivity between species or ecosystems. Ultimately, this dynamic system strives to establish equilibrium between tectonic uplift and river erosion. Determining whether or not a river network is in equilibrium, and, if not, what changes are required to bring it to equilibrium, will help us understand the processes underlying landscape evolution and the implications for river ecosystems.

### Methods

We developed the use of a proxy, referred to as χ, for steady-state river channel elevation. This proxy is based on the current geometry of the river network and provides a snapshot of the dynamic state of river basins. Geometric equilibrium in planform requires that a network map of χ exhibit equal values across all water divides (the ridges separating river basins). Disequilibrium river networks adjust their drainage area through divide migration (geometric change) or river capture (topologic change) until this condition is met. We constructed a numerical model to demonstrate that this is a fundamental characteristic of a stable river network. We applied this principle to natural landscapes using digital elevation models to calculate χ for three, very different, systems: the Loess Plateau in China, the eastern Central Range of Taiwan, and the southeastern United States.

### Results

The Loess Plateau is close to geometric equilibrium, with χ exhibiting nearly equal values across water divides. By contrast, the young and tectonically active Taiwan mountain belt is not in equilibrium, with numerous examples of actively migrating water divides and river network reorganization. The southeastern United States also appears to be far from equilibrium, with the Blue Ridge escarpment migrating to the northwest and the coastal plain rivers reorganizing in response to this change in boundary geometry. Major reorganization events, such as the capture of the headwaters of the Apalachicola River by the Savannah River, are readily identifiable in our maps.

### Discussion

Disequilibrium conditions in a river network imply greater variation of weathering, soil production, and erosion rates. Disequilibrium also implies more frequent river capture with implications for exchange of aquatic species and genetic diversification. Transient conditions in river basins are often interpreted in terms of tectonic perturbation, but our results show that river basin reorganization can occur even in tectonically quiescent regions such as the southeastern United States.

## River Reorganization

As rivers flow, they slowly transform the landscape. Their channels migrate, banks erode, and rivers can even flow backward if merged with another river. To understand how river basins balance erosional forces with regional tectonic uplift, **Willett et al.** (10.1126/science.1248765) analyzed maps of a proxy for river elevation and horizontal movement of river drainage divides across three large river systems in China, Taiwan, and the United States. Along with numerical modeling, the results demonstrate the degree to which these basins are at topographic equilibrium. The changing connectivity and topography of river networks influences how species migrate and how much material is delivered to larger bodies of water.

## Abstract

River networks evolve as migrating drainage divides reshape river basins and change network topology by capture of river channels. We demonstrate that a characteristic metric of river network geometry gauges the horizontal motion of drainage divides. Assessing this metric throughout a landscape maps the dynamic states of entire river networks, revealing diverse conditions: Drainage divides in the Loess Plateau of China appear stationary; the young topography of Taiwan has migrating divides driving adjustment of major basins; and rivers draining the ancient landscape of the southeastern United States are reorganizing in response to escarpment retreat and coastal advance. The ability to measure the dynamic reorganization of river basins presents opportunities to examine landscape-scale interactions among tectonics, erosion, and ecology.

The recognition that Earth’s surface topography evolves over time (*1*) provides a basis for reconstructing past tectonic (*2*, *3*) or climatic processes (*4*, *5*), as well as the formation of sedimentary sequences. The creation of geographic barriers or connections as topography evolves also influences speciation (*6*, *7*) and biodiversity (*8*). Despite the importance of landscape change and the long history of research, establishing rates of change, or even determining whether a given landscape is in a steady or disequilibrium state, has proven remarkably difficult. Geochronological techniques have made it possible to measure long-term rates of erosion (*9*, *10*), but large spatial and diverse temporal scales make it very challenging to assess topographic change by measuring erosion rates throughout an entire landscape.

The evolution of topography is fundamentally coupled to changes in river channel networks, which carve the landscape into ridges and valleys as they erode through rock. River basins are part of a dynamic system in which channel geometry, channel gradient, and network topology adjust toward a balance between tectonic uplift and erosion (*11*), a balance that depends on climatic conditions and the erodibility of the rock substrate (*12*). Examples of major changes in the topology or geometry of river networks have been reported on the basis of irregular network shape (*13*, *14*) or geologic evidence such as sediment provenance (*15*), and geologists have described these events with marked terms like river capture and stream piracy (*13*, *16*). However, without a straightforward way to map the motion of individual drainage divides, let alone the adjustment of an entire landscape, most studies of river systems have focused on the equilibration of vertical incision rates while assuming that map-view basin geometry remains fixed in time (*2*, *3*).

Recent efforts to generalize the long-term evolution of river networks have found support for drainage basin reorganization. Theoretical (*17*), experimental (*18*), and modeling studies (*19*–*22*) have argued that water divides are dynamic features of a landscape that routinely migrate, either progressively or through discrete river capture, in some cases even leading to complete reorganization of river networks (*17*). However, this theoretical progress has not been matched by observational support, which has remained limited to interpretations of large-scale capture events where geologic and geomorphic evidence is available (*14*, *23*, *24*).

We incorporate this view of mobile drainage divides into a procedure for mapping the dynamic evolution of river networks. We posit that equilibrium in a landscape requires that river network topology and drainage basin geometry adjust to maintain stationary water divides. This principle leads to a quantitative criterion for landscape disequilibrium that provides a prediction of drainage divide migration and a tool to map past, present, and future changes in a river network. We demonstrate this technique through analysis of a numerical model of landscape evolution and then apply it to several study sites, mapping transient features associated with river network reorganization.

## A Measure of River Basin Disequilibrium

Most landscapes comprise river channels that erode bedrock and transport sediment out of the system and hillslopes that erode by a variety of processes but remain in long-term equilibrium with their bounding river channels (*25*). The map-view arrangement of river basins evolves through the horizontal motion of drainage divides, which is a consequence of differential rates of river channel erosion on opposite sides of the divides. A measure of the disequilibrium of opposing river channels can therefore be used to infer the stability of the intervening divide and the direction of divide motion that would bring the adjoining channels toward equilibrium.

In a river network, two flow paths that originate at a common drainage divide and terminate at a common base level experience the same drop in elevation, but their steady-state elevation profiles—the theoretical profiles for which erosion would balance rock uplift—may differ depending on the topologic and geometric structure of the network. The uppermost parts of the flow paths, between the drainage divide and the channel heads where the river system initiates, traverse hillslopes. However, several factors imply that hillslopes do not contribute substantially to differences in the flow paths’ steady-state elevation profiles. First, the elevation drop along hillslopes is typically small relative to that in river channels. Second, steady-state hillslope relief is about equal on either side of the divide: Equilibrium hillslopes have similar gradients (*26*), and channel heads tend to be equidistant from divides (*25*). We therefore limit our analysis to the river channels. If drainage basins on opposite sides of a divide are structured such that the steady-state elevation of one channel head is higher than that of its opposing channel, the divide is not stable, and the network must adjust to bring the divide to a stable position. In general, any change in network topology or geometry that brings the cross-divide difference of steady-state channel head elevation to zero will bring the system into equilibrium. The simplest way to achieve this is by divide migration toward the basin with the higher steady-state elevation. Our strategy is to determine steady-state channel elevations and, thereby, directions and patterns of divide motion that allow the landscape to attain an equilibrium for both channel elevations and divide positions.

Steady-state channel elevation can be estimated from a model for river incision into bedrock. The most common model is based on the hypothesis that incision rate is proportional to the rate at which the flow expends energy as it travels over the bed, or stream power (*27*). Stream power depends on channel slope and river discharge, although the latter is generally replaced by the upstream drainage area, *A*, which serves as a readily available proxy. Including rock uplift due to tectonics, *U*, the elevation of a point in a channel, *z*, varies with time, *t*, and distance along a channel, *x*, according to (1)where *K* is a constant that depends on rock erodibility, precipitation rate (*28*), and channel geometry, and *m* and *n* are empirical constants. Other physical models, for example, based on bed shear stress, lead to an identical form of the erosion law; only the parameters, in particular, the exponents, are different. In most cases, these are determined empirically (see Materials and Methods), so Eq. 1 can be regarded as general for any erosion law that incorporates power law scaling between channel slope and drainage area at steady state. For the simple case of *U* and *K* constant in space and time, the steady-state solution of Eq. 1 is(2)where *z*_{b} is the elevation at the river network’s base level at *x* = *x*_{b}. The quantity χ is an integral function of position in the channel network (*29*),(3)where *A*_{0} is an arbitrary scaling area, and the integration is performed upstream from base level to location *x*. χ is the characteristic parameter for transient solutions of the linear (*n* = 1) version of Eq. 1 (*30*), and it remains the fundamental scaling parameter for the nonlinear case. The inclusion of the scaling area, *A*_{0}, gives χ dimensions of length, but the kinematic wave nature of Eq. 1 implies that χ could equally well represent a time. In particular, if *K ^{n}A*

_{0}

^{m}is included in the denominator of the integrand, χ takes on dimensions of time and, for the case of

*n*= 1, it becomes the characteristic time required for a perturbation at the river’s base level to reach a point

*x*in the channel (

*12*).

The term in parentheses in Eq. 2 represents the relative magnitudes of tectonic forcing and erosivity, and scales the magnitude of elevation. The parameter χ characterizes the river network topology and geometry, which determine how tectonic forcing generates variable topography throughout a river basin. Given the linear form of Eq. 2, it is apparent that χ serves as a metric for the steady-state elevation of a channel at location *x*. Thus, with constant tectonic forcing and homogeneous physical properties, a difference in χ across a divide implies disequilibrium and, presumably, motion of the divide in the direction of larger χ to achieve equilibrium (Fig. 1). This observation is the basis for our subsequent analysis: Mapping χ throughout a channel network and comparing χ values across drainage divides yield a snapshot of the dynamic reshaping of drainage basins.

## Elevation-χ Scaling with Changing Drainage Area

As a divide moves, either by continuous migration or through discrete river capture, drainage area is removed from one basin and added to the other. The channel length of each affected tributary also changes, leading to a change in the steady-state elevation of each channel head bounding the moving divide, presumably moving the channels toward equilibrium as in Fig. 1. However, analysis of a simplified scenario—the effect of a sudden change in drainage area on an equilibrium elevation profile (see Materials and Methods)—illustrates a feedback between erosion rate and divide motion that complicates this system. An instantaneous change in area induces an instantaneous change in χ, throwing the affected profile into a state of disequilibrium. Figure 2 shows the change in the χ plot (elevation against χ) of the perturbed channel for a given fractional increase or decrease of the upstream area. Area gain shifts the χ plot to the left, above the steady-state trend, and increases its length and thereby its maximum χ value, whereas area loss shifts the profile to the right, below the steady-state trend, and decreases its length and maximum χ value. A channel that lies above the steady-state trend on a χ plot erodes faster, on average, than the tectonic uplift rate (*29*), so a channel gaining area experiences an increase in average erosion rate, whereas a channel losing area experiences a decrease in average erosion rate. Branches of the channel network that connect to the affected channels do not necessarily experience any change in channel length, but they do experience the indirect effect of the change in erosion rate that propagates throughout the basin. This defines an important positive feedback in the system: A transfer of drainage area from one basin to another leads to changes throughout the affected drainage basins that encourage motion of the entire perimeters of the basins in the same direction as the original perturbation. The ultimate configuration of drainage divides if and when a landscape reaches equilibrium depends on the nonlinear interactions of multiple adjacent drainage basins and cannot easily be predicted. Here, we focus only on the local direction of divide motion toward equilibrium, but we identify some situations in which the positive feedback appears to dominate.

## Spatial Variations in Uplift Rate, Runoff, or Rock Erodibility

If *U* or *K* varies in space, and these variations are known, the solution for elevation can still be obtained by integration of Eq. 1. In practice, however, *U* and *K* are seldom known. It is more common to have information about relative values or spatial patterns. For example, uplift rate may vary across a fault; precipitation and runoff, which are included in *K*, may have a persistent spatial pattern; or rock erodibility may vary with rock type with a spatial distribution known from geologic mapping. If we express the spatial pattern of uplift and rock erodibility in terms of nondimensional functions of space, *U*^{*} and *K*^{*}, we can bring this variability inside the definition of χ without changing its dimensionality. Defining the uplift and erodibility as *U* = *U*_{0}*U*^{*}(*x*) and *K* = *K*_{0}*K*^{*}(*x*), where *U*_{0} and *K*_{0} are dimensional scaling values, the steady-state elevation is given by(4)with a modified χ defined as(5)

In landscapes with known spatial variations in *U* or *K*, Eq. 5 can be used to calculate χ′. However, if spatial variations in *U* and *K* are unknown but are likely to be minor, or if *U* and *K* vary in a systematic way over the river basins being analyzed, it may be possible to identify patterns of river basin reorganization using Eq. 3 to calculate χ, and interpreting spatial patterns of variation.

## Numerical Model of Landscape Evolution and χ Evolution

To test our hypothesis that mapping of χ characterizes the transient state of an evolving river network, we conducted an experiment using a numerical landscape evolution model (*21*). This model solves Eq. 1 for an evolving channel network spanning a grid of nodes and calculates the positions of water divides between channels using analytical solutions for the topography of low-order channels and hillslopes (see Materials and Methods).

There are many ways that river basin reorganization can be induced in a landscape. Here, we start with an initial steady topography generated using a linear gradient in the tectonic rock uplift rate from a lower value on the lower boundary to a higher value on the upper boundary (Fig. 3). This tectonic forcing leads to a highly asymmetrical mountain range in which the main drainage divide is closer to the upper boundary. To force river basin reorganization, we remove the uplift gradient, so that the entire domain experiences uniform uplift, and let the topography evolve until a new, symmetric steady-state mountain range develops.

The corresponding χ map (Fig. 3 and movie S1) confirms that the disequilibrium is characterized by differences in χ across divides, particularly the main divide. As predicted, divides move away from channels with low χ toward channels with high χ. In addition, we frequently observe isolated, low-order channels with high χ (Fig. 3), apparently in response to loss of drainage area and the positive feedback mechanism discussed above. These channels are transient features that tend to disappear quickly, suggesting that the area loss feedback eventually leads to a topological change in the network. Finally, we see that divides do not stabilize until χ values on either side are equal. At steady state, all points in the domain have the elevation predicted by χ, and there are no cross-divide changes in χ for adjoining channel heads. We also note that the time to steady state is many times longer than the profile equilibration time of the longest rivers in the model, demonstrating that the evolution of the network topology determines the equilibration time of the system.

As this example demonstrates, mapping χ across a landscape and noting discontinuities across drainage divides should reveal not only whether that landscape is in equilibrium but also the spatial pattern of divide migration, from which we can infer a history of landscape evolution and its forcings. We investigate this through a study of a series of natural examples below.

## Equilibrium Drainage Basins: China’s Loess Plateau

China’s Loess Plateau is a region in the Yellow River (Huang He) drainage basin where wind-blown sediments have accumulated to thicknesses of tens to hundreds of meters over the past 2.6 million years (*31*), draping preexisting topography (*32*). We constructed χ maps for two adjacent tributaries of the Yellow River, the Yanhe River and the Qingjian River (see Materials and Methods), which together drain an area of nearly 12,000 km^{2} (Fig. 4). There are no large contrasts in χ across drainage divides at any scale, with one exception in the headwaters of the Yanhe, where a strongly curved channel has lower χ than its neighbor to the north, a sign that it is gaining area by migration in this direction. Excepting this feature, the otherwise concordant χ values across drainage divides in this landscape suggest that the channel network is nearly stationary. This river channel network was established through the Tertiary on the tectonically stable Ordos platform, and the absence of tectonic perturbation has allowed it to attain a topologic (but not necessarily topographic) steady state. This may be additionally promoted by recent incision of the relatively uniform, erodible loess, although the main features of the drainage network appear to predate deposition of the loess (*32*).

## Basin Reorganization in an Active Mountain Range: Taiwan

The island of Taiwan was formed by collision of the Luzon volcanic arc with the Eurasian plate starting a few million years ago and continuing to the present day. Despite its youth, it has been argued that the overall height and width of the mountain range are in a steady state (*33*), in large part because the rates of tectonic shortening, uplift, and erosion are high (*34*). However, this condition does not necessarily require equilibrium at the scale of individual drainage basins. To investigate this point, we map χ in a limited area of the eastern Central Range (Fig. 5), where rock type, precipitation, and uplift rate are nearly uniform, and variations that do occur are primarily in an east-west direction, parallel to the dominant flow direction of the rivers. Rivers draining to the west of the main divide traverse different rock types and active structures with spatially variable uplift rates, so we avoid comparisons of χ across the main divide and limit our analysis to the internal and lateral divides of the eastward-draining rivers.

We find that many drainage divides at multiple scales have bounding channels with contrasting χ, which implies that the divides are unstable (Fig. 5). The high relief of Taiwan allows for easy visual comparison between cross-divide contrasts in χ and topographic features visible in topographic models and satellite imagery (Fig. 5, B and C, and fig. S4). In most cases, cross-divide contrasts in χ correspond to topographic features consistent with divide migration, including lower valley elevations, steeper channels, and evidence of faster hillslope erosion in basins predicted to be aggressors.

Disequilibrium is particularly evident for several long, isolated tributaries that have high χ relative to neighboring basins. Some of these high-χ tributaries have no side channels, indicating that the lateral divides are moving inward, starving side channels of drainage area, and transforming basin side slopes into unchanneled surfaces (for example, Fig. 5C). These tributaries appear to be victims of the positive feedback associated with area loss, in which loss of drainage area lowers erosion rate and increases χ, leading to surface uplift relative to neighboring basins, and increased vulnerability to further area loss.

There is also evidence for larger-scale basin reorganization. Figure 5B shows a basin that appears to be growing in all directions at the expense of its neighbors. The growing basin has steep upper channels and many landslide scars, indicating rapid channel incision. Both neighboring rivers appear to be losing drainage area in their upper reaches based on the existence of perched channels with high χ values. This suggests that even at the largest scale, the major river basins of the eastern Central Range appear to be adjusting their relative sizes, changing shape, and perhaps even increasing their number through expansion of smaller basins.

Comparison of χ and topography across divides also serves to confirm that there are no large changes in *K* or *U* through the study area. If this were the case, we would not expect to see the high degree of internal consistency between drainage divide asymmetry and cross-divide χ. Although we do recognize some instances where lithologic effects change channel steepness, these appear to be of minor importance compared to the geometric disequilibrium reflected in χ.

## Ongoing Basin Reorganization on a Passive Margin: Southeastern United States

The southeastern seaboard of the United States comprises the coastal plain and Piedmont physiographic provinces and is defined by the old mountain belt of the Appalachians in the west and by the Atlantic Ocean in the east, which formed by rifting starting about 200 million years ago (*35*). The boundary between the Piedmont and the Appalachians is the Blue Ridge (*36*, *37*), which defines the mountain front from Georgia to Pennsylvania and is interpreted to be an erosional escarpment associated with Atlantic rifting, similar to escarpments on other rifted margins (*38*). The top of the Blue Ridge defines the water divide between the Atlantic coastal rivers and the Mississippi River basin. This divide appears to be migrating inland, often by discrete capture of upland rivers (*15*, *23*, *36*, *39*). The spatial distribution of χ illuminates the escarpment beautifully (Fig. 6), with a strong discontinuity in χ between high values in the upper Tennessee River and Ohio River basins and low values in all basins draining directly to the east coast. This discontinuity defines an erosional front propagating to the northwest.

Some care must be taken in interpreting contrasts in χ across the main Appalachian divide, given that rivers that meet at this divide follow very different paths to the ocean and may therefore traverse regions with different uplift, climate, or rock type. If these variations were known precisely, they could be included in the characteristic parameter (Eq. 5). However, even without precise constraints on the spatial patterns of *U* and *K*, it is possible to show that the χ difference across the main divide is largely a signature of geometric disequilibrium. The most important heterogeneity that potentially creates spatially variable *K* is the contrast between the hard, erosion-resistant metamorphic and sedimentary rocks of the Appalachians and the softer sedimentary rocks of the Piedmont and coastal plain (*40*). If this effect were included through Eq. 5, the result would be to enhance the contrast in χ across the main divide, not diminish it.

Uplift in this tectonically quiescent region is dominated by isostatic response to erosion, which, if not constant, only varies over long wavelengths. However, there may be an additional, relatively young component of dynamic topography (*41*). This dynamic component may create a spatial gradient in *U*, but this gradient is oriented predominantly northeast-southwest along the Appalachian crest in the study area (*41*) and, hence, should not contribute to discontinuities in χ across the main divide. As with the Taiwan example, we can also compare more subtle features of the χ map and topography to check for consistency. For example, the divides between the Apalachicola River, the Tennessee River, and the Altamaha River (Fig. 6) show different values or even signs in cross-divide χ, and in all cases, the contrast is mimicked by subtle but clear asymmetry in topography. These observations suggest that the major features of the χ map do not arise from spatially variable *K* or *U*.

We also find support for inferred divide motion in erosion rates derived from cosmogenic nuclide concentrations in river sediment (*10*). Concentrations of ^{10}Be generated by cosmic ray exposure of quartz, now found as sand in the modern rivers, have been used to estimate erosion rates in a number of catchments across the southern Appalachians. We have compiled published data from the Great Smoky Mountains of Tennessee and North Carolina (*42*) and two segments of the Blue Ridge escarpment in North Carolina and Virginia (*43*) (Fig. 6). Motion of a drainage divide must ultimately be driven by different erosion rates on opposite sides of the divide. By comparing erosion rates from adjoining basins, we can estimate the direction of divide motion and test whether it is consistent with the direction inferred from the cross-divide difference in χ (see Materials and Methods). We found that the cross-divide χ difference correctly predicts the difference in erosion rate in 29 of the 34 adjoining basins we studied, with a strong correlation in magnitude in addition to predicting the correct sign (Fig. 6, inset).

Mapping of χ is also useful for identifying discrete capture events, many of which have been recognized along the Appalachian front. We present a documented example in Fig. 7 to show the signature pattern in χ. The Savannah River is proposed to have captured the headwaters of the Apalachicola River. Evidence for this interpretation includes the deeply incised Tallulah Gorge above the capture point and a common freshwater fauna in both catchments (*39*, *44*). The spatial distribution of χ shows sharp discontinuities across the Savannah-Apalachicola and Savannah-Tennessee divides, implying ongoing advance of the Savannah into both of these catchments (Fig. 7). In contrast, the Apalachicola-Tennessee divide is close to equilibrium. A discrete river capture involves a sudden transfer of drainage area from one basin to another, which will strongly affect χ but not elevation, thereby pushing the relationship between χ and elevation in opposite directions for the captured and capturing rivers (Figs. 1 and 2). This is demonstrated clearly in the Savannah capture region by tributaries to the upper Savannah, which gained area and therefore have smaller χ than expected for their elevation, and by the Apalachicola, which lost area and therefore has larger χ than expected for its elevation (Fig. 7). At present, the upper Savannah is responding to the increase in its drainage area by incising rapidly, creating marked features like the Tallulah Gorge as it lowers its elevation profile back to the regional χ-*z* trend. Similarly, the Apalachicola is incising less rapidly to raise its elevation profile back to the regional trend. This process will continue until both rivers reach equilibrium or, more likely, another capture occurs.

Another interesting feature of the map in Fig. 6 is the set of chevron-shaped drainage divides scattered across the coastal plain and Piedmont from Florida to Virginia, which are highlighted by consistently higher χ values on the coastal sides of the divides (Fig. 8). The anomalously high χ values in these basins are associated with anomalously high elevations (Fig. 8), similar to the area-starved basins in Taiwan (Fig. 5), but in this case, the topography has relief of only meters to tens of meters. The coastal basins are losing area as their headwaters narrow and shorten, and some may eventually disappear. This appears to be part of a regional basin reorganization in which the large Atlantic-draining basins are widening and pinching out the smaller, intervening basins. Such a process is consistent with northwest-southeast lengthening of the Piedmont and coastal plain provinces. River basins exhibit a narrow range of length-to-width ratios (*45*), so that lengthening in the flow direction can lead to basin reorganization, decreasing the number of basins to increase the average width (*17*, *18*, *46*). The Atlantic coastal region is lengthening by retreat of the Blue Ridge escarpment, as well as by seaward advance of the coastline, as evidenced by uplifted Pliocene terraces, possibly as a dynamic response to deep mantle flow (*41*).

## Global Implications

Our analysis presents and quantifies a dynamic view of landscapes. Drainage divides migrate or make discrete jumps; river basins expand, contract, and deform; and this dynamic reorganization can persist for hundreds of millions of years, or perhaps indefinitely in the presence of active tectonic deformation.

This evidence of shifting drainage divides helps explain well-known morphologic properties of drainage basins. Basins in varied geographic settings are observed to take on specific geometric forms or statistical states (*45*–*47*). In particular, characteristics such as length-area scaling, tributary junction angles, and fractal dimension appear to be globally consistent. Models that optimize local or global energy dissipation reproduce many of these fundamental network characteristics (*48*, *49*), although natural basins seem to be slightly suboptimal by these energy metrics. It has long been suspected that divide migration and river capture might be mechanisms by which drainage basins approach statistically uniform geometry (*50*). Our numerical models and analyses of natural landscapes provide evidence that this may, in fact, be the case, and additionally suggest that many natural basins are suboptimal because they have not yet reached an equilibrium configuration. Furthermore, our studies suggest that it is the geometric characteristics of the channel networks and their bounding divides that drive the system toward equilibrium.

The ability to map disequilibrium in river basins also has implications for other fields of study. For example, river profiles are often interpreted in terms of transient changes in tectonic uplift, climatic conditions, or sea level. However, our analysis suggests that many perturbations in river profiles may instead arise from changes in drainage area. Drainage network connections and water divides that form geographic barriers affect the transport of sediment, nutrients, and dissolved elements, including those with important global biogeochemical functions such as nitrogen and carbon. Similarly, the migration or diversification of aquatic species and entire ecosystems depends on the transport pathways defined by rivers. Our ability to identify past and ongoing changes in river networks creates a new opportunity to explore connections between geological, chemical, and biological systems.

## Materials and Methods

### Response of χ to a Change in Drainage Area

For the calculation in Fig. 2, we assume that area and channel length scale according to *A* = *k*_{a}*x ^{h}* with

*x*= 0 at the water divide (

*51*), and to keep the calculation analytical, we assume that

*h*= 2 and

*m*/

*n*= 0.5. Most other values for the coefficients

*h*,

*m*, and

*n*would require numerical integration, but the results will not change in character for a reasonable range of these coefficients. Defining a nondimensional length , we obtain an initial distribution of(6)and a steady-state elevation linearly proportional to χ

_{init}. We consider an instantaneous change in area of Δ

*A*of this basin, assumed to be the consequence of adding or removing area to the headwaters of the basin. The change in area can be positive or negative, but it must occur upstream of the analysis. The width of the basin downstream of the area change is unchanged. With the nondimensional area perturbation, , following the area perturbation is defined by(7)which can be integrated to give(8)

Plotting χ_{init} against χ_{p} (Fig. 2) is equivalent to showing elevation against χ for an instantaneous change in area. For a basin losing area, there is a corresponding decrease in length of the channel, hence the smaller extent in χ_{p} of the curves in Fig. 2. For the basin gaining area, there is a corresponding increase in channel length, but this is not shown in Fig. 2 because the actual slope and length of this segment depend on the area distribution of the captured channel. A plot of channel slope against drainage area shows a similar perturbation (fig. S1).

### Landscape Evolution Model

The numerical model for landscape evolution (Fig. 3) uses the code DAC described by Goren *et al.* (*21*). DAC incorporates numerical and analytical solutions to represent processes at different scales. The numerical component solves the conservation of mass equation with a stream power incision law (Eq. 1) for a channel network spanning a dynamic, irregular grid. An analytical solution is used for the subgrid topography to represent low-order fluvial channels and hillslopes. The channel network and the numerical grid evolve as divide migration leads to river capture, abandonment of channels with insufficient drainage area, or creation of new channels as hillslopes lengthen. The model in Fig. 3 consists of a 50 × 100–km rectangular domain in which the four edges are fixed to a constant elevation. Precipitation and rock type are steady and uniform. An initial topography is generated by imposing a rock uplift rate that varies linearly from 0.5 mm/year along the lower edge to 5.0 mm/year along the upper edge. These conditions are run to steady state. Subsequently, the tectonic gradient is removed, and uplift rate is set to a constant value of 1 mm/year. Basin divides on all scales become unstable and migrate, causing drainage basins to change their size and shape (see movie S1). χ is calculated from Eq. 3 on the discrete grid, but because part of the river channel network is contained in the analytical solution, sections of the uppermost catchments are not shown in Fig. 3 or movie S1.

### Construction of χ Maps

Construction of maps of χ followed a specific protocol. Hydraulic attributes of base level, flow directions, flow paths, and accumulated flow (upstream drainage area) were extracted from a digital elevation model (DEM). For flow direction and paths, closed basins were filled, and steepest descent neighbors were found for local flow direction. Any pixel with a contributing area less than a critical value, *A*_{c}, was excluded from the analysis. The critical area was typically on the order of 10^{6} m^{2}, although it can vary depending on the DEM resolution and the drainage density of the landscape. This value does not affect the downstream value of χ, but it does determine how high into a catchment we conduct the analysis. An arbitrary scaling area, *A*_{0}, and a value for *m*/*n* were selected, followed by integration of Eq. 3 to determine χ for all pixels in the domain.

The concavity, *m*/*n*, was selected through an iterative process. First, we constructed a series of χ plots for individual drainage basins, covering a range of concavity values, and noted the *m*/*n* value that best reduced the scatter (*29*, *52*). This was used to construct an initial χ map, which was then interpreted in terms of divide stability. If divides appeared stable within the basin of interest, the process was complete. However, if there were large contrasts in cross-divide χ, we reconsidered the χ plots, noting the profiles of channels that appeared to be gaining or losing area based on the χ analysis, or on independent geologic, geomorphic, or geochemical data. Because these χ plots are expected to curve up or down (Figs. 1 and 2), we modified our selection of *m*/*n* accordingly. The interpretation of area gain or loss is partially based on the resultant χ map, so is potentially circular, but the iterative process consistently converged to a single solution. In some cases, we also applied this method to individual channel profiles where we could make an a priori assumption regarding the shape of the χ plot. As a final check, we overlaid the χ map on the digital topography and checked that predicted cross-divide χ discontinuities were consistent with topographic features such as asymmetric upper channel profiles across divides. Figures 5 and 8 are examples of such comparison.

### Loess Plateau

The 90-m resolution DEM derived from the Shuttle Radar Topography Mission (SRTM) (*53*) was used for the Loess Plateau. Two drainage basins were extracted and analyzed using a critical area, *A*_{c}, of 0.05 km^{2} and a scaling area, *A*_{0}, of 1 m^{2}. The low value for the critical area reflects the high drainage density of the area. Interpretation of the resultant χ-elevation plots gave an optimum *m*/*n* of 0.35 (fig. S2).

### Taiwan

For Taiwan, we used a 40-m resolution DEM derived from aerial photographs and available from the Center for Space and Remote Sensing Research, National Central University, Taiwan. A scaling area, *A*_{0}, of 1 m^{2} and a critical area, *A*_{c}, of 1 km^{2} were used. Plots of χ-elevation were constructed with concavity varying from 0.0 to 0.6. The smallest scatter is obtained for a concavity of 0.3 to 0.4, but we favor a slightly higher value. Given our finding of many moving divides and captures, we expect many channels to exhibit χ plots with low slope and high χ, as well as many channels with kinks in the χ-elevation profile in response to capture. These features become more obvious with concavities of 0.45 to 0.55 (fig. S3), so we prefer this value and use a value of 0.5 for the maps of χ. This value was confirmed by comparison with topography and imagery, which is quite marked in Taiwan (Fig. 3 and fig. S4).

### Southeastern United States

For the analysis of the southeastern United States, we used the 90-m CGIAR SRTM DEM, analyzed with a critical area, *A*_{c}, of 0.3 km^{2} and a scaling area, *A*_{0}, of 1 m^{2}. We constructed χ-elevation plots for individual drainage basins (identified in fig. S5) in the study area for a range of *m*/*n* values (figs. S6 to S9). We found very dynamic divides in the region with examples of basins that were dominantly growing and others that were dominantly shrinking. The minimum variance in the χ-elevation plots was found for values of *m*/*n* from 0.25 to 0.35, but to match the χ-elevation plot concavity to our interpretation of growing and shrinking basins, we required a larger *m*/*n* value of 0.4 to 0.5. We also conducted a similar analysis on individual rivers where geological or geomorphic evidence permitted an a priori assessment of capture or a moving divide. This exercise also gave *m*/*n* values of 0.4 to 0.5 (for example, fig. S10). Maps in this paper were constructed using an *m*/*n* of 0.45, though other values were constructed to test sensitivity to *m*/*n* (fig. S11), which revealed that the features described in this paper were robust over a large range of *m*/*n* values.

Map construction on the western side of the Appalachians is complicated by the need to include the drainage area of the entire Mississippi River basin. To avoid integrating the alluvial lower Mississippi, we initiated the χ integration at the confluences between the westward draining rivers and the Ohio River. This required selecting a χ-elevation pair for each river. To calculate these values, we used the regional χ-elevation trend for the lower Tennessee River and lower Kanawha River. This resulted in an initial point of 8 m in χ and 100 m in elevation for the Tennessee River and 10 m in χ and 148 m in elevation for the Kanawha River.

### Differential Erosion Rates Estimated from Cosmogenic Radionuclide Concentrations

We used published concentrations of the CRN ^{10}Be in quartz river sand to estimate differential erosion rates in select river basins. We used published data from Portenga and Bierman (*10*) in three areas within the southern Appalachians: two locations along the Blue Ridge escarpment, and one in the Great Smoky Mountains (Fig. 6). Original studies are in (*42*, *43*). We selected locations with CRN measurements in adjoining basins where the difference in erosion rate can be compared to the cross-divide contrast in χ. To estimate differential erosion across a divide, we differenced the basin-averaged erosion rates of adjoining basins and assigned the resulting differential erosion rate, Δ*E*, to the divide segment common to both basins. We then calculated the mean Δχ across the shared divide segment by differencing pairs of adjoining channel heads in the same direction. This procedure was carried out for all primary divide segments between basins with comparable size. We did not analyze nested sub-basins. We also excluded one basin from the Great Smoky Mountains location, which appeared to have anomalous ^{10}Be concentrations associated with an unusual intrusive igneous unit (figs. S12 to S14 and table S1).

In the two Blue Ridge locations, there were few basins with CRN data directly adjoining one another. In these locations, we calculated an average differential erosion rate across the escarpment itself. We divided the escarpment into three segments (two segments as shown in fig. S14, and one segment shown in fig. S13) that had erosion rate measurements on each side. We then determined the mean erosion rate of each side of the escarpment by taking the mean erosion rate of basins proximal to the escarpment, but draining in opposing directions. We differenced these mean erosion rates to calculate an escarpment-averaged Δ*E*. We then calculated the mean Δχ across each escarpment segment by differencing, in the same direction, the mean channel head χ values along the escarpment segment. The results from all three study areas were normalized by dividing Δ*E* and Δχ by their respective mean *E* or χ for the two basins being compared (Fig. 6, inset).

We also tested how much of the overall variance in the CRN erosion rate data could be explained by the cross-divide χ gradients. We calculated an average cross-divide Δ_{χ} for a basin by integrating the difference of channel head χ across the exterior perimeter of each basin. This provides an aggressivity metric, in which positive values indicate basins that are likely growing at the expense of their neighbors, whereas basins with negative values are victims that are likely to be shrinking. These aggressivity indices are compared with basin-wide erosion rates in fig. S15.

**Data Archive.** Chi maps for each region can be downloaded in kml format from www.sciencemag.org/content/343/6175/1248765/suppl/DC1. Data details are given in table S2.

## Supplementary Materials

www.sciencemag.org/content/343/6175/1248765/suppl/DC1

Figs. S1 to S15

Tables S1 and S2

Movie S1

Databases S1 to S9

## References and Notes

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
**Acknowledgments:**S.D.W., S.W.M., and J.T.P. thank the Canadian Institute for Advanced Research (CIFAR) for research support and members of the Earth System Evolution Program of CIFAR for discussions. S.W.M. was supported by a CIFAR postdoctoral fellowship and by the U.S. National Science Foundation through award EAR-0951672 to J.T.P. C.-Y.C. was supported by the Swiss National Research Foundation. Data are archived in the supplementary materials.