Watersheds in disordered media

What is the best way to divide a rugged landscape? Since ancient times, watersheds separating adjacent water systems that flow, for example, toward different seas, have been used to delimit boundaries. Interestingly, serious and even tense border disputes between countries have relied on the subtle geometrical properties of these tortuous lines. For instance, slight and even anthropogenic modifications of landscapes can produce large changes in a watershed, and the effects can be highly nonlocal. Although the watershed concept arises naturally in geomorphology, where it plays a fundamental role in water management, landslide, and flood prevention, it also has important applications in seemingly unrelated fields such as image processing and medicine. Despite the far-reaching consequences of the scaling properties on watershed-related hydrological and political issues, it was only recently that a more profound and revealing connection has been disclosed between the concept of watershed and statistical physics of disordered systems. This review initially surveys the origin and definition of a watershed line in a geomorphological framework to subsequently introduce its basic geometrical and physical properties. Results on statistical properties of watersheds obtained from artificial model landscapes generated with long-range correlations are presented and shown to be in good qualitative and quantitative agreement with real landscapes.


I. INTRODUCTION
Although both start in the mountains of Switzerland, the Rhine and Rhone rivers diverge while flowing towards different seas. While the Rhine empties into the North Sea, the Rhone drains into the Mediterranean. Similarly, all over the world one finds rivers with close sources but distant mouths. For example, the Colorado and the Rio Grande even open towards different oceans. When looking at a landscape, how to identify the regions draining towards one side or the other? When rain falls or snow melts on a landscape, the dynamics of the surface water is determined by the topography of the landscape. Water flows downhill and overpasses small mounds by forming lakes which eventually overspill. A drainage basin is then defined as the region where water flows towards the same outlet. Its shape and extension strongly depend on the topography. The line separating two adjacent basins is the watershed line, which typically wanders along the mountain crests [1].
Since watersheds provide information about the dynamics of surface water, they play a fundamental role in water management [2][3][4], landslides [5][6][7][8], and flood prevention [8][9][10]. Watersheds are also of relevance in the political context they have been used to demarcate borders between countries such as the one between Argentina and Chile [11]. Thus, the understanding of their statistical properties and resilience to changes in the to- * nmaraujo@fc.ul.pt † kjs73@cam.ac.uk ‡ hans@ifb.baug.ethz.ch § soares@fisica.ufc.br pography of the landscape are two important questions that we will review here. In general, for every landscape, several outlets can be defined, each one with a corresponding drainage basin. Sets of small drainage basins eventually drain towards the same outlet forming a even larger basin. Such hierarchy results in a larger number of watersheds. However, without loss of generality, the study of watershed lines can be simplified by splitting the landscape into only two large basins, each one draining towards opposite boundaries. Figure 1 shows how this watershed can be identified by flooding the landscape from the valleys (see supplemental video). Two sinks are initially defined (the lower-left and upper-right boundaries in the example). While flooding the landscape, each time two lakes (A and B) draining towards opposite sinks are about to connect, one imposes a physical barrier between the two. In the end, the watershed line is the line formed by the barriers that separate these two lakes.
The concept of watersheds is also of interest in other fields like, e.g., in medical image processing [12]. There, computed tomography scans need to be segmented to identify different tissues. The pictures are discretized into pixels and a number is assigned to each one of them according to the intensity. The segmentation procedure consists in clustering neighboring pixels following the order of increasing intensity gradient, splitting the image into different parts (tissues). The equivalent to the watershed line corresponds to the line separating two different tissues [13,14].
It was recently shown that watersheds can be described in the context of percolation theory in terms of bridges and cutting bond models [15], with numerical evidence that they are Schramm-Loewer Evolution (SLE) curves in the continuum limit [16]. This association explains why watersheds on uncorrelated landscapes as well as other statistical physics models, such as, optimal path cracks [17][18][19], fuse networks [20], and loopless percolation [15], belong to the same universality class of optimal paths in strongly disordered media.

II. THE LANDSCAPE
The motion of surface water is defined by the topology of the landscape, usually characterized by the spatial distribution of heights. Although this distribution is a continuous function for real landscapes, it is typically coarse-grained and represented as a digital elevation map (DEM) of regular cells (e.g., a square lattice of sites or bonds) to which average heights can be associated. This process is exemplarily shown in Figs. 2(a)-(c). In fact, modern procedures of analyzing real landscapes numerically process Grayscale Digital Images where the gray intensity of each pixel is transformed into a height, resulting in a DEM. As such, discretized maps have been useful to delimit spatial boundaries in a wide range of problems, from tracing watersheds and river networks in landscapes [21][22][23][24][25] to the identification of cancerous cells in human tissues [13,26], and the study of spatial competition in multispecies ecosystems [27,28].
To study watersheds theoretically one can generate artificial landscapes. Starting with a regular lattice, a numerical value is randomly assigned to each element, cor- responding to its average height. If these values are spatially correlated, a correlated artificial landscape is obtained. Otherwise, the landscape is called an uncorrelated artificial landscape. Natural landscapes are characterized by long-range correlations [29]. In Sec. VII we will review how to generate correlated artificial landscapes and their main properties. The DEM can be further simplified by mapping them onto ranked surfaces, where every element (site or bond) has a unique rank associated with its corresponding value [15]. A ranked surface can be defined in the following way. Given a two-dimensional discretized map of size L × L, one generates a list containing the heights of its elements (sites or bonds) in crescent order, and then replaces the numerical values in the original map by their corresponding (integer) ranks. As depicted in Fig. 2(d), the result is a ranked surface.

IV. FRACTAL DIMENSION
Breyer and Snow [35] have studied 12 basins in the United States and concluded that their watershed lines are self-similar objects [36,37]. A self-similar structure is characterized by its fractal dimension d f , which is defined here through the scaling of the number M of edges belonging to the watershed with the linear size L of the discretized map, They obtained fractal dimensions in the range 1.05−1.12.
Fehr and co-workers [32], with the method described in Sec. III confirmed this self-similar behavior over more than three orders of magnitude and measured the fractal dimension for watersheds in several real landscapes, as summarized in Tab. I. For uncorrelated artificial landscapes the watershed fractal dimension has been estimated to be d f = 1.2168 ± 0.0005 [32,38,39]. This value has drawn considerable attention since it was also found in several other physical models such as optimal paths in strong disorder and optimal path cracks [17-19, 40, 41], bridge percolation [15,38,42], and the surface of explosive percolation clusters [43,44]. The conjecture that all these models might belong to the same universality class has opened a broad range of possible implications and ap-plications of the properties of watersheds. As discussed in Ref. [15], the relation between most of such models can be established when they are described within the framework of ranked surfaces (see also Sec. II).

V. WATERSHEDS IN THREE AND HIGHER DIMENSIONS
Up to now, we have focused on the watershed line that divides the landscape into drainage basins for water on the surface. However, in reality, water also penetrates the soil and flows underground. Thus, the concept of watersheds and ranked surfaces can be extended to three dimensions. The soil can be described as a porous medium consisting of a network of pores connected through channels. When a fluid penetrates through this medium, a threshold pressure p k can be defined for each channel k such that the channel can only be invaded when p ≥ p k , where p is the fluid pressure. In general, a channel k is closed when p k > p, and open otherwise. This system can be mapped into a three dimensional ranked volume, where the lattice elements are the channels and the rank is defined by the increasing order of the threshold pressure. The sequence in the rank corresponds to the order of channel openings when the fluid pressure is quasistatically raised from zero. Analogously to the ranked surfaces, one can split the space in two regions, draining towards opposite boundaries.
The watershed in three dimensions is a surface of fractal dimension d f = 2.487 ± 0.003 [39]. An example for a simple-cubic ranked volume is shown in Fig. 5. In general, for lattices of size L d , where d is the spatial dimension, the watershed blocks connectivity from one side to the other. Thus, the watershed fractal dimension must follow d − 1 ≤ d f ≤ d, i.e., with increasing dimension d f also increases [15].

VI. IMPACT OF PERTURBATION ON WATERSHEDS
The stability of watersheds is also a subject of interest. For example, changes in the watershed might affect the sediment supply of rivers [45]. Also, the understanding of the temporal evolution of drainage networks provides valuable insight into the biodiversity between basins [46]. Geographers and geomorphologists have found that the evolution of watersheds is typically driven by local changes of the landscape [47]. These events can be triggered by various mechanisms like erosion [46,48,49], natural damming [8], tectonic motion [50][51][52], as well as volcanic activity [53]. Although rare, these local events can have a huge impact on the hydrological system [47][48][49]54]. For example, it was shown that a local height change of less than two meters at a location close to the Kashabowie Provincial Park, some kilometers North of the US-Canadian border, can trigger a displacement in the watershed such that the area enclosed by the original and the new watersheds is about 3730 km 2 [55].
The stability of watersheds is also relevant in the political context. For example, Chile and Argentina share a common border with more than 5000 kilometers, which was the source of a long dispute between these countries [11]. A treaty established this border as being the watershed between the Atlantic and Pacific Oceans for several segments. In 1902, the Argentinian Francisco Moreno contributed significantly to elucidate the technical basis for dispute. He proved that during the quarternary glaciations, the watershed line changed. In particular, several Patagonian lakes currently draining to the Pacific Ocean were in fact originally part of the Atlantic Ocean basin. Consequently, he argued, instead of belonging to Chile they should be awarded to Argentina.
Douglas and Schmeeckle [45] have performed fifteen table top experiments to study the mechanisms of drainage rearrangements. In spite of being diverse, the mechanisms triggering the evolution of watersheds are all modifications of the topography [8,46,50]. In that perspective, the effect of such events can be investigated by applying small local perturbations to natural and artificial landscapes and analyzing the changes in the water-shed [34,55]. Specifically, one starts with a discretized landscape and computes the original watershed. A local event is then induced by changing the height of a site k h k → h k + ∆, where ∆ is the perturbation strength, and the new watershed is identified. The impact of perturbations on watersheds in two and three dimensions will be discussed in the following.

A. Impact of perturbations in two dimensions
Fehr and co-workers [55] have normalized the perturbation strength by the height difference between the highest and lowest height of the landscape. For each landscape, they have sequentially perturbed every site of the corresponding discretized map, such that each perturbed landscape differs from the original only in one single site. The effect of those perturbations affecting the watershed has been quantified by the properties of the region enclosed by the original and the new watersheds. For that region, they have measured its area, corresponding to the number of sites, N s , in the discretized map and the distance R between its original and new outlets, defined as the points where water escapes from this region. Note that, there are only two outlets involved in this procedure: one related to the original landscape and the other to the new one. Since ∆ is strictly positive, the outlet in the new landscape corresponds always to the perturbed site k.
Scale-free behavior has been found for the distribution P (N s ) of the number of enclosed sites N s , the probability distribution P (R) of the distance R between outlets, and the dependence of the average N s on R [34,55]. Specifically, where the measured exponents are summarized in Tab. II. The power-law decay (2b) and the relation (2c) imply that a localized perturbation can have a large impact on the shape of the watershed even at very large distances, hence having a non-local effect. Additionally, the analysis of the fraction of perturbed sites affecting the watershed revealed a power-law scaling with the strength ∆. This finding supports the conclusion that changes in the watershed can be even triggered by anthropological small perturbations [55]. For the region enclosed by the original and the new watersheds, an invasion percolation (IP) cluster can be obtained by imposing a pressure drop between the outlet in the new watershed and the one in the original one, always invading along the steepest descent of the entire cluster perimeter. The size distribution P (M IP |R) of these clusters, for each fixed distance R between outlets, has been shown to scale as, where α * ≈ 1.39. This exponent corresponds to the one found for the size distribution of point-to-point IPclusters [57]. In that process, invasion clusters are obtained in a random medium by invading from one point to another at a certain distance. By contrast, in the watershed case the invasion is always started from the outlet on the new watershed. This difference between starting at any point or in the outlet justifies the additional factor of unity in the scaling (Eq. (3)), since the probabilities need to be rescaled by the size of the IP-cluster [55].

B. Impact of perturbations in three dimensions
The impact of perturbations on watersheds has been also analyzed in three dimensions [34]. Similar to 2D, a perturbation is induced by changing the value of a site k, h k → h k + ∆, where ∆ is the perturbation strength. The number of sites N s enclosed by the original and new watersheds corresponds to a volume. Power-law scaling in terms of Eqs. (2a)-(2c) has also been observed with the exponents summarized in Tab. II. The value of α * is similar to the one found in two dimensions. According to Lee [58], the size distribution of the point-to-point IP-cluster is independent on the dimensionality of the system. Therefore, the numerical agreement between α * for different spatial dimensions supports the relation with invasion percolation.
Fehr et al. [34,55] used fractional Brownian motion (fBm) on a square lattice [71] to incorporate long-range correlations controlled by the Hurst exponent H. They calculated how the fractal dimension decreases with the Hurst exponent and found good quantitative agreement with the exponents obtained for natural landscapes, typically with 0.3 < H < 0.5, which is the known range of Hurst exponents for real landscapes on length scales larger than 1 km (see Ref. [29] and references therein). They also obtained α, β and ρ for several values of H, observing that both β and ρ increase with H, finding also good quantitative agreement. Thus, their model provides a complete quantitative description of the effects observed on natural landscapes.

VIII. FINAL REMARKS
Here we solely discussed cases with one watershed, but the same theoretical framework can be straightforwardly extended to tackle other space partition problems. For example, the identification of the entire set of watersheds of a landscape with multiple outlets helps identifying the catchment areas contributing to each river or reservoir [72]. A systematic study of disordered media with multiple outlets is still missing. Examples of open questions are: How does the distribution of catchment basins or the number of triplets (points where two watersheds meet) depend on the number of outlets? And, how does the statistics of perturbations change in the presence of triplets? The division of a volume into several parts is also a problem of practical interest in the extraction of resources from porous soils [73]. Studies of such systems have mainly considered uncorrelated disordered media. The role of long-range correlation there is still an open problem.