While the hierarchical organization of river basins has been well described, little physical evidence exists on how such signatures change in time as a function of drainage network evolution. Using a soil-mantled experimental landscape, a fourth-order rill network was developed in response to applied rainfall and base-level lowering. Drainage basin area and tributary length scales were determined at discrete times during the evolution of the rill network. Results show that rill basin area exhibited self-similarity across rill order (space) and as the network grew in time, and that the coefficient of Hack’s law, once the network began to grow, also was invariant with time. This experimental evidence strongly supports the use of scaling arguments for drainage basin evolution and for the prediction of rill network development and organization on hillslopes and agricultural fields.


The areal dimensions of river basins and the lengths of the longest tributaries have been shown to exhibit scaling relationships (Rodríguez-Iturbe and Rinaldo, 1997; Rinaldo et al., 2014). Drainage basin area, A, can be rendered as the planar projection of its width, L, and length, L||, or A = sLL|| where s is a basin shape factor (Rigon et al., 1996). If the ratio L/L|| and s remain constant for all A, then basin area is considered self-similar. If L/L|| increases or decreases but s remains constant with increasing A, then basin area is considered self-affine. The tendency for river basins to become longer and narrower with increasing area, and to exhibit self-affinity, has been suggested by Hack’s law (Hack, 1957), such that LAh where L is the length of the longest stream and h ≈ 0.6. Evidence also suggests that the magnitude of Hack’s exponent could vary as a function of drainage-basin scale and that this variation reflects the dominant geomorphic processes in operation (Dodds and Rothman, 2000; Gangodagamage et al., 2011). Abundant field data support the tendency for river basins to approximate these scaling relationships (Rigon et al., 1996; Dodds and Rothman, 2000), and these observations have informed theory and numerical models for landscape evolution driven by fluvial erosion (Rodríguez-Iturbe et al., 1992; Veneziano and Niemann, 2000; Pelletier, 2007).

Yet little physical evidence currently exists on an evolving drainage system and how these hierarchical signatures change in time as a function of drainage network development. That is, these scaling relationships have been derived for static river landscapes rather than dynamic or evolving landscapes. While never stated explicitly, space or scale may have been used as a surrogate for time; small basins within a watershed could represent earlier stages of basin evolution. Moreover, soil erosion prediction technology used on hillslopes and agricultural fields excludes the contribution of individual branches of evolving rill networks in sediment discharge determinations (Govers et al., 2007). While models are available to predict flow and sediment transport within individual rills (Papanicolaou et al., 2010), no technology is currently available to predict the development, extension, and hierarchical organization of an evolving rill network.

Physical models of rill networks, however, can be used to address these important gaps in the time evolution of self-formed drainage systems (Schumm et al., 1987; Pelletier, 2003; Bennett et al., 2015a). Our paper explores the evolution of rill networks on an experimental landscape exogenically forced by rainfall, overland flow, and base-level adjustments. Details of the experiments, and the morphodynamics of the evolving rill networks, are presented elsewhere (Gordon et al., 2012; Bennett et al., 2015b). The focus here is on defining tributary basin area scales and Hack’s law and on tracking their temporal variation as a function of rill network development. It will be shown that the incised rill tributary basins exhibit self-similarity relatively early in the development of the drainage system, and that the rill basins maintain these scaling relationships during their subsequent growth and evolution.


The experiments were conducted using a flume 7 m long, 2.4 m wide, and 0.3 m deep, which had diversions at the downstream end to funnel water and sediment through a 0.16-m-wide outlet (Gordon et al., 2012; Bennett et al., 2015b). Bed slope was set at 5%, and the perforated flume bottom was covered with a 0.02 m layer of 10 mm gravel and landscape fabric to allow for free drainage through the soil bed. A sandy loam (20% clay, 22% silt, and 58% sand by mass) was air dried, passed through a 2 mm sieve, placed in the flume in 0.02 m layers, and packed to a target bulk density of ∼1360 kg/m. The final soil bed ranged in thickness from 0.15 m in the center to 0.25 m along the sidewalls, and the bed surface was fashioned to assume a given initial topography (Fig. 1A). At the basin outlet, a collection trough was attached to an adjustable broad-crested weir that allowed the flume base level to be dropped incrementally.

Simulated rainfall at an average rate of 87 ± 31 mm/hr was administered to the soil surface under pressure by cone-jet nozzles at a fall height of 3.6 m, and 20 storms were applied to the prepared bed. Following the initial 85 min storm, the base level at the downstream weir was lowered 0.03 m, which created a pre-formed step at the flume outlet that proceeded to migrate through the landscape once rainfall application and surface runoff recommenced. Subsequent storms had varying durations depending upon how much rill erosion occurred. Base level was lowered a second time after 850 min, which produced another wave of rill incision. Flow and sediment exiting the flume were collected in 1 L bottles, the water was decanted, and the sediment was oven dried and weighed to determine rates of sediment discharge QS (kg/s) and runoff Q (m3/s; Fig. 2).

Digital elevation models (DEMs) were produced using photogrammetric techniques. Images of the bed surface were obtained using a 10 megapixel digital camera attached to a moveable carriage suspended above the flume. Commercially available software (Leica Photogrammetry Suite v. 9.2) was used to create DEMs at a user-specified resolution of 3 mm (Figs. 1A and 1B), with horizontal and vertical uncertainties on the order of 1 and 3 mm, respectively (see Chandler et al., 2001; Gordon et al., 2012). GIS techniques were employed to delineate and define drainage and rill network characteristics, accomplished using ArcGIS 10.1 (http://www.arcgis.com; with an accumulation threshold of 5000 grid cells); these results were in agreement with those of TauDEM (http://hydrology.usu.edu/taudem/taudem5/index.html; Tarboton, 1997) and GeoNet (Passalacqua et al., 2010).


Response of Soil-Mantled Surface to Exogenic Forcing

An incised rill network developed relatively quickly in time in response to rainfall application and base-level lowering. Runoff from the soil surface began after 18 min and rose to an asymptotic rate of ∼2.69 × 10−4 m3/s within 80 min of runtime (Fig. 2). This initial runoff produced minor erosion and sediment discharge, most notably in the central portion of the landscape (Fig. 1B). Each base-level adjustment produced a wave of degradation and landscape incision manifested by the upstream migration of a primary headcut, whose scour depth approximated the downstream adjustment and which bifurcated into newly created lower-order tributaries, and by spikes in sediment discharge that decreased asymptotically with time (Fig. 2; Gordon et al., 2012). As a result of this exogenic forcing, a fourth-order rill network developed (Fig. 1B). The time variation in drainage density, defined as Dd = ∑ L/A where L is the Euclidian distance of all rills, shows an asymptotic growth when total flume area AF is used and an asymptotic decline when the incised rill network area AR is used (Fig. 2).

Basin Area Elasticity, Hack’s Law, and Rill Evolution in Time and Space

As the landscape evolved, the characteristics of the rill tributary basins were determined for discrete times (100, 120, 180, 230, 370, 850, and 1060 min after the beginning of rainfall application) where tributary orders were identified and drainage basin areas were demarcated. The planar projection of length L|| and width L of each basin, as defined by the basin drainage area, then was objectively determined for all rill orders using the minimum boundary geometry tool (rectangle-by-width method) in ArcGIS. An example of this procedure is shown for the 1060 min surface (Fig. 1C). The results of this automated method were compared directly to boundary boxes drawn by hand (the 18 second-order [2°] basins in Fig. 1C), and the scale ratios L/L|| for these two procedures were statistically indistinguishable (p < 0.001, using a t-statistic for the Pearson’s correlation coefficient; Rogerson, 2006).

Basin boundary boxes were defined for all rill tributaries for the times chosen, and the ratios of width to length L/L|| for each order and for each time, as well as metrics of the distribution, are summarized in Figure 3. As the rill network evolved from the 100 min to the 1060 min time steps, the number of tributary rills increased (except for the 4° rill): from 18 to 69 for 1°, from 4 to 18 for 2°, and from 1 to 6 for 3°. The ranges in L/L|| provide evidence for the inherent elasticity in the basin length scales for a rill network freely developed on an erodible soil. Yet the values of L/L|| show clear central tendencies. The mean L/L|| values for all 1°, 2°, and 3° rills range from 0.335 to 0.567 (Fig. 3), with an ensemble mean value of ∼0.480 ± 0.171. This is in general agreement with observations made in rivers (Rigon et al., 1996).

Statistical tests confirm that the average basin width-to-length ratio does not change with rill order at a given time or in time for a given rill order. Excluding all data with less than four observations, the distributions of L/L|| were analyzed statistically to determine (1) if the mean values varied as a function rill order for a given time (i.e., 1° versus 2° at 120 min, 2° versus 3° at 120 min, etc.; performed using the nonparametric Kruskall-Wallis test; Rogerson 2006), and (2) if the mean values varied as a function time for a given order (i.e., 1° at 100 min versus 1° at 120 min, 1° at 100 min versus 1° at 180 min, etc.). For the comparison of L/L|| as a function of rill order for a given time, no mean values were statistically different (p >> 0.05 for all tests). For the comparison of L/L|| as a function of time for a given rill order, no mean values were statistically different (p >> 0.05 for all tests). By definition, using L/L||, all 1°, 2°, and 3° rill tributary basin areas demonstrate self-similarity (1) in space for a given time during network evolution, and (2) in time as the rill network evolved in space. To assess the reproducibility of these results, the same network analysis described above was performed on a companion experiment (run 1 at 787 min; Gordon et al., 2012). The mean value of L/L|| for all rill basins in this experiment was close to 0.5, and the mean values of L/L|| as a function of rill order (1°, 2°, and 3° rills) for this given time were not statistically different, corroborating the results herein.

The coefficient h in Hack’s law also was derived for each time increment of the evolving rill network (Fig. 2). Values range from 0.288 to 0.533, with an average of 0.461 ± 0.069. Statistical tests performed on these results confirm that h does not vary with time once the rill network begins to develop. That is, h = 0.288 ± 0.076 at 100 min, and this value is statistically different from all subsequent coefficients derived for the evolving rill network (p < 0.025 for most comparisons; t-statistic for regression coefficients; Rogerson, 2006). But all values of h after the 100 min time interval are statistically indistinguishable (p >> 0.05), with an average value of 0.491 ± 0.067, similar in magnitude to L/L|| noted above.


For a soil-mantled experimental landscape exogenically forced by simulated rainfall and base-level lowering, the resultant rill network displays self-similarity in planar basin area across rill order or space (for the freely formed 1°, 2°, and 3° rills) and with time, with L/L|| ≈ 0.48. Application of Hack’s law also shows that h ≈ 0.49 soon after initial network development, which is consistent with the basin area self-similarity observation. For the single 4° rill whose drainage area effectively includes the entire landscape, L/L|| ≈ 0.371 (Fig. 3), which is close in magnitude to the that of the experimental facility (0.343). This observation suggests that the scale of the largest network could be affected by the boundary conditions imposed. These experimental observations represent the first physical evidence for scaling ratios of an evolving drainage network, in which drainage area increased by a factor of 16, drainage density increased by a factor of six, and the number of tributary links increased by a factor of four.

The broadest implications of this work lie in the prediction of drainage network development and landscape evolution of hillslopes. Self-similarity of rill network evolution in these ideal experimental conditions provides strong support for the lateral and headward growth (expansion) of a space-filling network unencumbered by material characteristics, stratification, or hydrologic heterogeneities (Howard, 1971; Schumm et al., 1987; Hancock and Willgoose, 2001). While the tributary basins display some elasticity in size, the self-similarity of basin area and the invariance of Hack’s law strongly suggest an intrinsic threshold for rill tributary length L|| for a given width L, so that L|| ≈ 2L everywhere within the network. That is, fingertip tributaries extending across a landscape will be limited by the lateral space available to maintain basin area self-similarity, potentially forcing the rill tributary to bifurcate during network growth and expansion. Such empirical evidence provides strong justification for the use of scaling arguments as a framework for the prediction of drainage basin evolution, and supporting theory already exists. Smith (2010) presented a model for the growth and evolution of river basins driven by surface runoff. This analysis showed that during the initial stages of landscape evolution, as the drainage networks are just emerging, the largest and fastest growing basins have L/L|| ≈ 0.5, consistent with the observations reported here. Lastly, soil erosion prediction technology at the plot scale continues to treat rill erosion as a lumped process (Govers et al., 2007) as opposed to discrete morphodynamic channels of concentrated flow (Papanicolaou et al., 2010). The empirical evidence presented here also provides strong support for the use of hierarchical scaling arguments necessary to predict the growth, evolution, and organization of rill networks on hillslopes and agriculture fields.


Using a soil-mantled experimental landscape subjected to simulated rain and base-level lowering, basin area and stream length scales were determined for an evolving rill network. The ratio of basin width to length was close to 0.5 for all rill orders and for all stages of network development. The coefficient for Hack’s law quickly achieved a value close to 0.5, and it too did not vary as the network continued to expand. These results demonstrate that as a rill network grows in size under optimum experimental conditions, rill tributaries display self-similarity in basin area in space (across rill orders) and in time. These observations provide strong support for the use of scaling relationships for the prediction of rill network evolution in time and space and their potential use in assessment of soil erosion and landscape degradation.

This research was supported by the National Science Foundation (grant BCS 1359904). Constructive comments on the manuscript offered by T. Smith and three anonymous referees are greatly appreciated. The authors would be happy to provide electronic copies of these experimental data to interested parties.