## Abstract

River network scaling laws describe how their shape varies with their size. However, the regional variation of this size-dependence remains poorly understood. Here we show that river network scaling laws vary systematically with the climatic aridity index. We find that arid basins do not change their proportions with size, while humid basins do. To explore why, we study an aspect ratio *L*_{⊥}/*L*_{∥} between basin width *L*_{⊥} and basin length *L*_{∥}. We find that the aspect ratio exhibits a dependence on climate and argue that this can be understood as a structural consequence of the confluence angle. We then find that, in humid basins, the aspect ratio decreases with basin size, which we attribute to a common hydrogeological hierarchy. Our results offer an explanation of the variability in network scaling exponents and suggest that the absence of self-similarity in humid basins can be understood as a morphological expression of subsurface processes.

## 1. Introduction

River networks form stunning branching patterns that repeat over many scales. In recent years, their seemingly fractal appearance has led to study in the context of simple unifying dynamics that operate independent of scale [1–3]. Perhaps the conceptually simplest of these is the random walk model introduced by Scheidegger [4], who suggested river flow direction could be approximated as a random downhill choice on a lattice. A number of other simple models have also been devised. For example, a simple loosening of the downhill condition of the Scheidegger model (a random walk on a lattice in one of eight directions, rather than two or three) produces a network with different scaling properties [5]. Similarly, a class of optimal channel network models [1,6] suggests that river networks self-organize to an energetically optimal state, bearing qualitative similarity to avalanches and other so-called self-organized critical processes [7]. Ultimately, study of such models revealed how scaling laws delineate distinct morphological classes of river networks [1,3,8,9]. However, real river dynamics and forms are often subject to a wide range of geological and environmental constraints [10] so that each river network exhibits different scaling laws according to its environment. These constraints inevitably introduce variability into network scaling laws and have stymied efforts to develop a unifying theory [1,3,11,12].

Here we ask how exponents of these network scaling laws can be understood through both an environmental gradient and the simple physical models described above. We then consider how these scaling laws can be understood explicitly through study of basin shape and its dependence on local aridity and a critical length scale associated with the presence of shallow groundwater.

## 2. Scaling, universality and climate

River basin scaling is frequently characterized through a canonical scaling law known as Hack's law [13], which relates the length of the longest stream in a network *l* to basin area *a* by an exponent *h*, known as the Hack exponent:
*h*. Namely, we calculate the aridity index over the contiguous United States, defined as *P* is precipitation rate and PET is potential evapotranspiration. High values of A correspond to humid regions, while low values of A correspond to arid regions. This calculation is described further in appendix A. We fit lines to data of log_{10}*l* against log_{10}*a* to obtain *h* for arid and humid basins, shown in figure 1. We find that arid basins have a Hack exponent of *h*≃0.5, while humid basins have *h*≃0.6, suggesting that the aridity index parameterizes the scaling exponent *h*.

Although Hack's law is canonically used to study how basin shape changes with scale, it alone does not provide a complete picture of river network scaling [3,16]. We must also consider the sinuosity of streams [12], which can be expressed through the relation:
*d* is a scaling exponent and *L*_{∥} is the basin length. Together, *d* and *h* define a *universality class* of river networks [3,8,9,16], and adherence to a particular universality class would suggest that a river network's attributes can be largely understood through an underlying mathematical model.

The behaviours of *d* and *h* with A are shown in figure 2*a,b*. Figure 2*c* compares *d*, *h*)≃(1.04, 0.50) in arid environments and (*d*, *h*)≃(1.09, 0.59) in humid environments. These exponents suggest that, relative to humid basins, river networks in arid environments are roughly self-similar, in that arid basins do not change shape with changing scale. Arid networks also resemble the ground state of one of the optimal channel networks studied in [9], whereas those in humid environments do not. Still, neither the arid nor humid basins bear a definite resemblance to these universality classes of river networks. We note that the values of *h* and *d* corresponding to optimal channel networks in figure 2 are exact results deriving from global minimization of an energy functional. However, numerical simulation of optimal channel networks typically result in local minimization, yielding *h* = 0.58 ± 0.03 [1].

## 3. Basin shape

The differences in the Hack exponent *h* for humid and arid basins should be manifested in the size-dependence of basin shape. We, therefore, explicitly consider basin shape to shed light on factors that may explain the scaling behaviour we have observed. We define a basin width *L*_{⊥} = *a*/*L*_{∥}. The aspect ratio *L*_{⊥}/*L*_{∥} then characterizes the shape of a basin. Figure 3*a* shows a map of *L*_{⊥}/*L*_{∥}. See appendix A for a full explanation of this calculation. We find that the aspect ratio *L*_{⊥}/*L*_{∥} correlates with aridity, as shown in figure 3*b*. In particular, basins in humid regions appear to be proportionally fatter than basins in arid regions.

To understand why, we consider how confluence angle can influence the shape of isometric basins. To this end, consider a network that can be structurally approximated as a confluence branching symmetrically with an angle *α* between confluent branches off of a longer parent branch, as illustrated in figure 4. We can naturally draw a boundary between two branches, and thus between the two adjacent basins, that is delineated by the direction of the parent channel. If we assume that these basins are symmetric (the area to the left of the stream is equal to the area on the right), then it follows that the drainage area *a* of a basin of longitudinal extent *L*_{∥} is simply a slice of angle *α* (radians) of a circle with radius *L*_{∥}:
*a*≃*L*_{⊥}*L*_{∥}, the aspect ratio

Recent work has demonstrated that groundwater-fed channel networks grow in a direction that is determined by the shape of the surrounding groundwater table, predicting a branching angle of 2*π*/5 = 72° [17,18]. And an analysis of the contiguous United States finds that this angle appears widespread in the local branching angles of humid environments [14], suggesting that the dynamics of groundwater processes may be influential in determining junction angle. Arid environments, however, exhibit a narrower branching angle approximately 45°. If these two end-member branching angles are related to network-wide structures, we can predict aspect ratios for the model network of figure 4 from the observed branching angles associated with humid and arid environments: a branching angle of 72° corresponds to *L*_{⊥}/*L*_{∥} = *π*/5 = 0.63, and a branching angle of 45° corresponds to *L*_{⊥}/*L*_{∥} = *π*/8 = 0.39. These values are indicated by solid blue and yellow lines in figure 3*b*, respectively. The rough agreement suggests that the branching angles can naturally produce the observed changes in basin shape. This supports the idea that groundwater plays a substantial role in shaping river basins in humid environments.

Finally, to return to the question of variability in the Hack exponent, we consider how basin shape changes with scale. Figure 5 shows *L*_{⊥}/*L*_{∥} as a function of basin area *a* for two different aridity ranges. We find that the aspect ratio in arid regions remains ≃0.4 for basins of all sizes, indicating that in arid basins, proportions do not change with scale. However, the aspect ratio in humid regions indeed decreases as basin size increases, suggesting that their proportions do change with scale, in agreement with our measurements of the basin scaling laws in figures 1 and 2. Specifically, humid basins appear to exhibit an aspect ratio of about *π*/5 = 0.63 when small, but as basin size increases, there is an apparent transition towards a smaller aspect ratio of about *π*/8 = 0.39 characteristic of arid basins.

## 4. Deeper connections

We hypothesize that this scale-dependent humid basin shape may be an expression of the scale-dependence of subsurface water flow. In particular, the groundwater processes that give rise to the 72° confluence angle [17] and the *π*/5 aspect ratio in humid basins require nonlocal communication between branches through the groundwater field, but this interaction may be limited in its spatial extent. To understand the limiting length of subsurface interaction, we require a characterization of the groundwater table and its relationship with the topography. To this end, we consider two lengths: the typical distance *L*_{s} between streams and the typical distance *L*_{w} between streams that intersect the water table. An illustration of these length scales can be found in figure 6. By definition, *L*_{s} ≤ *L*_{w}. When *L*_{s}∼*L*_{w}, most streams present are groundwater-fed. When *L*_{w}≫*L*_{s}, a number of streams form from overland flow and are decoupled from subsurface flow.

Inspired by references [19,20], we define *L*_{w} as
*R* is the recharge rate, *k* is the hydraulic conductivity, *H* is the thickness of the water table, *D* measures topographic relief and *m* is a constant. A derivation motivating this quantity can be found in appendix B.

We calculate *L*_{w} throughout the contiguous United States using compilations of hydraulic conductivity, recharge rate and elevation following ref. [20] (see appendix A). Figure 7*a* shows the geographic distribution of *L*_{w}, and figure 7*b* shows *L*_{⊥}/*L*_{∥} as a function of basin size for *L*_{w} < 10 m (blue), 10 m ≤*L*_{w}≤100 m (green) and *L*_{w} > 100 m (yellow). Although we do not know the detailed geographic variation of *L*_{s}, it likely ranges between 0.1 and 1 km [21,22]. Therefore, *L*_{w} < 100 m (blue) likely corresponds to *L*_{w}∼*L*_{s} and *L*_{w} > 100 m (yellow) likely corresponds to *L*_{w} > *L*_{s}. In the latter case, the water table does not intersect low-order streams. The aspect ratio of basins beyond this scale decreases to less than about 0.5, consistent with figure 5.

On the other hand, when *L*_{w} < 100 m, *L*_{w}∼*L*_{s} and streams are fed by groundwater. In this case, we find that the aspect ratio of small basins is ∼*π*/5, but it becomes smaller with increasing basin size. Its behaviour with scale can be understood through the flow regimes described by Toth [23]. On the local (and smallest) scale, the shape of the groundwater table is not strongly correlated to small fluctuations in the local topographic surface, suggesting that the growth of streams on this scale is dominated by local groundwater flow [24]. It is on this scale that planform geometry can be understood through the groundwater field [17,25,26]. However, on the larger, regional scale, the water table replicates the topography and the resulting flows should follow topographic gradients [19,27]. Distant streams are still coupled by groundwater flow, but the transit times of these groundwater flows are long and the fluxes are small [28,29]. The decay from the wide aspect ratio in figure 7*b* is therefore gradual, not immediate.

The convergence of the humid basin aspect ratio towards the arid aspect ratio may indicate that large, humid basins self-organize in the same way as arid basins of all sizes. This convergence is consistent with earlier observations that large basins are self-similar, independent of climate and geology [30,31].

## 5. Summary

We have shown that arid basins have a Hack exponent of 0.5, suggestive of a self-similar model of growth in which network proportions do not change with scale. However, humid basins do not exhibit such self-similarity, growing thinner as they grow larger. We argue that this size-dependent humid basin shape results from a physical limit to groundwater processes, and that the observed Hack exponent of *h*≃0.6 is a manifestation of this limit. Remarkably, these relations are found despite the natural variability one expects in a variety of geologic settings.

Our work suggests that humid river basins exhibit a unique, mathematical dichotomy. At the scale of large basins, arid and humid networks behave alike, collecting water, perhaps doing so according to an optimality principle [1,32]. Yet, on the smallest scale, humid basins grow in response to a diffusive groundwater field [17,26]. In this way, the extent to which the water table is submerged below the surface determines the geometry of drainage networks at the surface.

## Data accessibility

Raw river data were obtained from NHDPlus v. 2, available at http://www.horizon-systems.com/NHDPlus/NHDPlusV2_home.php.

## Author's contributions

R.S.Y., E.S. and D.H.R. designed the research; R.S.Y., A.A.P., E.S. and D.H.R. performed the research; R.S.Y., A.A.P. and H.S. analyzed data; R.S.Y. and D.H.R. wrote the paper.

## Competing interests

We declare we have no competing interests.

## Funding

This material is based upon work supported by the US Department of Energy, Office of Science, Office of Basic Energy Sciences, Chemical Sciences, Geosciences and Biosciences Division under award number FG02-99ER15004.

## Acknowledgments

We thank Yossi Cohen and Olivier Devauchelle for their insightful comments.

## Appendix A

We study the geometry of drainage basins across the contiguous United States, derived from the NHD Plus v. 2 data set [33]. Each drainage basin is comprised of a number of sub-basins, each with a total upstream drainage area *a* and longest stream length *l*. We find the longitudinal basin length, *L*_{∥} by taking the distance between the start and end coordinates of the longest stream. These start and end coordinates were obtained by projecting the NHDPlus stream networks onto the Lambert conformal cone 102004 projection as described in [14]. To obtain the aspect ratio *L*_{⊥}/*L*_{∥}, we observe that because *a*≃*L*_{⊥}*L*_{∥}, we can write this aspect ratio equivalently as *a*/*L*^{2}_{∥}, which we then calculate for all basins and sub-basins within the contiguous United States.

We calculate the Hack exponent *h* in figure 2 by obtaining the slope of a line fit by least squares to log_{10}*l* against log_{10}*a*. We similarly calculate the sinuosity exponent *d* by fitting log_{10}*l* against log_{10}*L*_{∥}. We produce the map of *h* by applying the same methods to all basins geographically contained in each HUC6 basin. For all other maps, we average values within each HUC6 basin. To determine the climatic dependence of these exponents, we calculate the aridity index *P* is the precipitation rate, and PET is the potential evapotranspiration. Within each sub-basin, we take the geometric mean of A over the entire length of all rivers in the sub-basin. For the calculation of the network scaling exponents, we used geographically diverse points which exhibited values of A within the plotted bin ranges. We calculate *L*_{w} using data from references [20,35,36],^{1} Gleeson *et al.* [20] and equation 4.1. We consider the topographic relief *D* to be the maximum terrain rise between neighbouring river junctions.

## Appendix B. The length scale *L*_{w}

Consider a region of terrain such that within the domain, the dominant form of water transportation is the transport of groundwater, but any groundwater that reaches the edge of the domain is drained away; call the length-scale of the domain *L*, such that the area of the domain is like *L*^{2} times a small constant. A simple example of such a domain might be a circular region of radius *L* with no rivers on the interior but rivers along most of the outside edge.

Suppose that the height of the groundwater at each point in the domain is *H* + *h*, where *H* is a large constant, and *h* accounts for the (presumably small) variations of height within the domain. Then by Darcy's law, the groundwater discharge at a point is
*k* is the hydraulic conductivity, and ∇*h* is the slope in *h*.

By reckoning the balance between water flowing in and out of the domain, we can estimate a typical magnitude for *h* within the domain necessary to make the flow of water balance. The total rate at which water is added to the groundwater in the domain through rainfall is equal to the *groundwater recharge rate* *R* times the area of the domain. This equals the total outflow of water along the border of the domain:
*R* is roughly constant over the whole domain and that the magnitude of ∇*h* is also roughly constant along the boundary, and we ignore constants to let *L*^{2} stand for the area of the domain and *L* to stand for its perimeter, we get
*h*| is *RL*/*kH*, times small constants. Thus the typical variation in groundwater height *h* over the lengthscale *L* of the domain is
*L* surrounded by water on all sides, with a recharge rate *R* that is constant within the whole domain, we can exactly solve the Poisson equation
*h* (in the middle) and the lowest value of *h* (on the border) is exactly
*L*_{w} = *L* and *D* = Δ*h* and recover equation 4.1 for the case *m* = 1 after rearrangement of equation B 4.

## Footnotes

↵1 For Tropical Agriculture IC Resampled SRTM v. 4.1. data retrieved from CIAT, see http://gisweb.ciat.cgiar.org/TRMM/SRTM_Resampled_250m/.

- Received February 6, 2018.
- Accepted June 13, 2018.

- © 2018 The Author(s)

Published by the Royal Society. All rights reserved.