Non-linear hydrologic organization

We revisit three variants of the well-known Stommel diagrams that have been used to summarize knowledge of characteristic scales in time and space of some important hydrologic phenomena and modified these diagrams focusing on spatiotemporal scaling analyses of the underlying hydrologic processes. In the present paper we focus on soil formation, vegetation growth, and drainage network organization. We use existing scaling relationships for vegetation growth and soil formation, both of which refer to the same fundamental length and timescales defining flow rates at the pore scale but different powers of the power law relating time and space. The principle of a hierarchical organization of optimal subsurface flow paths could underlie both root lateral spread (RLS) of vegetation and drainage basin organization. To assess the applicability of scaling, and to extend the Stommel diagrams, data for soil depth, vegetation root lateral spread, and drainage basin length have been accessed. The new data considered here include timescales out to 150 Myr that correspond to depths of up to 240 m and horizontal length scales up to 6400 km and probe the limits of drainage basin development in time, depth, and horizontal extent.

1 Introduction

The development of “physically based” and verifiable spatiotemporal scaling relationships is one of the important goals stated in the publication Opportunities in the Hydrologic Sciences (National Research Council, 1991), which helped guide the establishment of NSF's Hydrologic Sciences Program in the Earth Sciences (EAR) Division. Our purpose is to find techniques which may help to develop and test such scaling relationships.

Hydrogeologists and earth scientists occasionally organize (eco)hydrologic phenomena according to their spatial and temporal scales so as to locate them in a single figure (National Research Council, 1991, 2001; Bloeschl and Sivapalan, 1995; Loague and Corwin, 2006). Such figures, known as Stommel (1963) diagrams, are illustrative and can, under the right circumstances, trigger further useful work, such as testing hypotheses regarding appropriate spatiotemporal scaling relationships and the relevant space and timescales that enter in. Further, the same figures may even be able to clarify basic tenets about surface compared with subsurface hydrology and link processes not originally known to be connected, even those that cross hydrologic boundaries.

In this particular note, we consider a few basic concepts presented in such figures in the light of recent work on the spatial scaling of river networks (Hunt, 2016, 2017a), the spatiotemporal scaling of chemical weathering (Hunt et al., 2014a; Hunt and Ghanbarian, 2016), soil formation (Yu and Hunt, 2017a, b; Yu et al., 2017; Egli et al., 2018), and vegetation growth (Hunt, 2017b). Guided by Stommel diagrams, we use these latter relationships to expand our spatial scaling of river networks to a spatiotemporal scaling framework.

Chemical weathering, which is typically limited by downward water flux (infiltration) in the unsaturated zone, is the principle limiting factor of soil development (Egli et al., 2018; Hunt et al., 2021). Thus, soil depth will have a dependence on that downward flux. Vegetation growth is related to transpiration (Hunt et al., 2020a) (Fig. 1). In each case, the fundamental network length scale is on the order of a particle (or plant xylem) size. The magnitude of groundwater fluxes, which are often nearly horizontal, is typically in the tens of meters per year (Bloeschl and Sivapalan, 1995), but magnitude of typical flux through the unsaturated zone and groundwater recharge is typically (e.g., Lvovitch, 1973) on the order of 10 %–30 % of the total precipitation (ranging from 0.001 to 10 m yr⁻¹). Infiltration and transpiration velocities are typically on the order from 0.1 to 1 m yr⁻¹, but extreme climate cases (e.g., the Atacama Desert and the New Zealand Alps) may extend this variability by orders of magnitude.

Figure 1 Growth rates and transpiration rates of eucalypts (i.e., data from Roberts et al., 2001; Givnish et al., 2014). The unknown parameter $x_{\mathrm{0}}/t_{\mathrm{0}}$ is estimated from the precipitation, which is variable in the study region but typically about 2 m yr⁻¹ where the trees were found. The 3D value of the optimal path exponent (shown) is in much closer accord (than the 2D value) with the data over the entire age range from 1 month to 300 years. The extracted parameters for transpiration and tree height as functions of time were not distinguishable (a ca. 4 % difference in individual parameter values).

Our proposed power-law length–time scaling relationship can be represented in the following form (Hunt, 2017b):

$$\begin{array}{}\text{(1)}& x=x_{\mathrm{0}}{\left({\displaystyle \frac{t}{t_{\mathrm{0}}}}\right)}^s.\end{array}$$

In this relationship, x is a distance (for soil the depth; for plant growth it is the root lateral spread, RLS; and for river networks it is the length of a river), t is the time, x₀ is the network characteristic scale (a pore separation), and the ratio of $x_{\mathrm{0}}/t_{\mathrm{0}}$ is the appropriate pore-scale flow rate, v₀. The power s was chosen (Hunt, 2017b), for soil development, as the inverse of the percolation backbone dimensionality, $\mathrm{1}/D_{\mathrm{b}}=\mathrm{0.53}$, and, for vegetation growth, equal to the inverse of the 2D percolation optimal path exponent, $\mathrm{1}/D_{\text{opt}}=\mathrm{0.83}$. These non-linear exponents reflect diminishing connectivity of the network with increasing length scales. For heterogeneous networks, D_opt gives the fractal dimensionality of the path connecting two (parallel) planar surfaces with the lowest possible cumulative resistance (Porto et al., 1997), while D_b is the fractal dimensionality of the percolation backbone, which is known to govern the spatiotemporal scaling of solute transport (Sahimi, 2014; Lee et al., 1999). Solute transport is, in turn, known to limit chemical weathering and soil formation (Yu and Hunt, 2017a, b; Hunt et al., 2021).

It was suggested that the variability in flow rates would be responsible for the most important variability in actual scaling relationships, though the fundamental length scale could also be of significance. The flow rate variability was assumed (Hunt, 2017b) to be ca. 0.024 to 20 m yr⁻¹, with the upper limit maybe too high, but which appears to give a good accounting for all three phenomena: soil formation, vegetation growth, and river drainage development. The relevance of Eq. (1) using the same parameters x₀ and t₀ from Yu et al. (2017a, b) and Hunt (2017b) will be tested further below by investigating its ability to frame the understanding of the Stommel diagrams.

2 Spatial and temporal scales of processes of relevance to hydrology

2.1 The Loague and Corwin (2006) Stommel diagram

In the context of our work on spatiotemporal scaling, we modified the Stommel diagram plotted in the chapter “Scale Issues” of the Handbook of Groundwater Engineering by Loague and Corwin (2006), which is shown in Fig. 2. We have added the colored strips to indicate important spatiotemporal scaling relationships of Eq. (1). The width of the strips indicates roughly 2 orders of magnitude variability in mean annual flow rates, which is largely a result of water supply and its variability from meteorological driving forces but also the heterogeneity of soil properties. We also added the icon for mineral dissolution and changed the position of landscape evolution so that it would represent a weathered depth.

Figure 2 Organization of hydrologic processes according to spatial and temporal scales with logarithmic axes (modified after Loague and Corwin, 2006). Thus, a constant slope on this figure represents the value of a power in a power law. The blue strip represents water flow distances for scale-invariant flow rates as a function of time and has a slope of 1. The green strip denotes root radial extent (and vegetation height) as a function of time and has a slope of 0.83. The brown strip represents soil depth as a function of time and has a slope of 0.53. These slopes less than 1 imply that vegetation physical extent and soil depth have decreasing growth rates with time. The particular exponents that govern these scaling relationships are the optimal path exponent in two dimensions for vegetation growth and the percolation backbone in three dimensions for soil formation. See the text for further detail.

The blue strip represents instantaneous pore-scale subsurface flow rates within an order of magnitude of 1 µm s⁻¹ (about 30 m yr⁻¹) (reported in a Stommel diagram of Bloeschl and Sivapalan, 1995, as the range of subsurface flow rates most commonly observed). This value is generally compatible with an estimate (US Geological Survey, 2021) that regional (mainly horizontal) groundwater flow rates range between 1 ft yr⁻¹ and 1 ft d⁻¹ (or between 0.3 and 100 m yr⁻¹). On a bilogarithmic plot, a process with constant velocity x=v₀t shows up with only a single relevant parameter (the velocity, v₀), which determines the vertical position; the range of likely velocities then generates the range of vertical positions. A flow rate as large as 1 µm s⁻¹ can be observed in saturated, more nearly horizontal flow, below the water table. However, the largest vertical flow rates in the unsaturated zone are likely to be as much as a factor 3 smaller. Nevertheless, a maximum water flow rate of ca. $\mathrm{0.63}\mathrm{\mu }\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}=\mathrm{20}$ m yr⁻¹ has proved (Hunt, 2017b) an excellent predictor of the maximum soil depth over time periods ranging from decades to greater than 10 Myr and is reused here.

In Fig. 2, we give the water flow rates as scale independent even though subsurface flow rates are assumed in Bloeschl and Sivapalan (1995) and inferred from observation (Skoien et al., 2003) to increase with scale. The question is of considerable importance, though discussions have not provided unambiguous results. This assumption already underlay the construction of the original figure in Hunt (2017b) and is given further basis in Hunt et al. (2014b).

Multiplication of the velocity of reaction products by their molar density yields their rate of reaction, when such reactions are transport-limited. This scaling velocity diminishes as a negative power of the time, according to the time derivative of Eq. (1) (Ghanbarian-Alavijeh et al., 2012; Hunt et al., 2014a; Hunt and Ghanbarian, 2016). Proportionality of the reaction rate to the soil formation rate allows for a prediction of the soil depth using Eq. (1) when erosion rates are negligible (for times less than required to reach steady state; Yu and Hunt, 2017a). Thus, the simple scaling relationship that defines the solute transport distance (Eq. 1) for advective solute transport in three-dimensional heterogeneous porous media (Ghanbarian and Hunt, 2012) together with the same range of flow rates yields the brown strip. In this Stommel figure, we have used the value of the instantaneous flow rate compatibly with both Loague and Corwin (2006) and Bloeschl and Sivapalan (1995). But, in order to make individual predictions of soil depth as a function of time (Yu and Hunt, 2017a, b; Egli et al., 2018; Yu et al., 2019), the flow rate for each field site where soil depths are measured must be chosen as the annual mean value for the vertical flow in the unsaturated zone, expressed as the quotient of P−ET and the porosity (P= precipitation, ET = evapotranspiration), which is typically an order of magnitude (or more) smaller than 20 m yr⁻¹. The vertical placement of the brown strip is chosen to make the ranges of solute and flow velocities equal (brown and blue strip widths equal) at the length scale corresponding to a typical silt-sized grain of about 30 µm (and timescale of around 30 s). Thus, the reduction of solute velocity with increasing scale sets on at the fundamental network scale of a single grain. Note that this band matches well with the icons of “chemical weathering”, “soil depth”, and “landscape evolution”, as well as reasonably well with “soil structure” and “mineral dissolution”. This suggests the interpretation, confirmed in a number of references (Yu and Hunt, 2017a, b; Yu et al., 2017, 2019; Egli et al., 2018), that use of the solute velocity in advective flow to predict chemical weathering and soil formation rates is accurate on the average to within 20 % (see, particularly, Egli et al., 2018, for verifications).

The green strip in Fig. 2 is defined through the optimal path exponent of percolation theory in two dimensions and represents equally total transpiration as well as RLS as a function of time (Hunt et al., 2020a). Optimal paths are defined by the minimal flow resistance. For periods up to a century, RLSs are nearly the same as tree height (Hunt and Manzoni, 2015, based on data from Phillips et al., 2014, 2015; Kalliokoski et al., 2008, and Stone and Kalisz, 1991). The vertical position of this band in Fig. 2 is also chosen so that the variability of the root tip extension rate at the timescale of about 30 s is equal to that of the water flow rate (blue strip) in the medium. The accompanying interpretation is that roots growing through heterogeneous soils tend to follow paths with lowest cumulative resistance. The access to nutrients near the surface (delivering the 2D topology assumed) provided by this growth pattern is essential for plant growth. Note that the green band appears to be consistent with the icons “grass” and “leaf canopy”, while it is also in accord with “vegetated buffer strip”. It is reasonable that “soil moisture” extends across the flow and vegetation growth bands, since the former position on the diagram represents the rate at which moisture is delivered to the soil from precipitation and the latter represents the dominant rate at which it is taken up by transpiration. The mathematical predictions for vegetation growth have been verified separately (Hunt, 2017b; Hunt et al., 2020a, b).

While a number of icons fall into the positions expected if such phenomena are indeed governed by the processes suggested, there are a few that fall notably outside the predicted ranges. “Crops” fall on (or below) the blue strip. The relatively minor limitation on crop growth to a large flow rate is provided by ample watering and fertilization, meaning that root growth in a nutrient/water search is not a limiting factor; indeed crop growth is mostly linear in time (Hunt, 2017b). “Forest” falls also below the blue strip but for a different reason. Reflection on the possible relevance of a length scale of 100 m or more in a year indicates that this length does not refer to RLS or the height of individual trees but more likely to the horizontal spread of forests into newly available habitats. This spread can be mediated by atmospheric processes, such as seed dispersal and, thus, is not subsurface limited. As expected, “weathering rind” thicknesses are less than “soil depths”, since, e.g., subsurface clasts used for measuring weathering rind thicknesses have much lower hydraulic conductivity values than the soil (Sak, 2021), and the amount of water entering the weathering rind is reduced by approximately the contrast in hydraulic conductivity values. Finally, the position of “soil texture” at around 30 µm after about 1000 years represents a typical grain size in the silt range, compatible with the fundamental scale of the network, though it could also be located at scales of centimeters, which defines a typical length scale on which heterogeneity of soil texture is defined.

An important principle revealed in Fig. 2 is its organization of surface and subsurface phenomena. The fundamental subsurface flow rate is not assumed to change with scale. Note that all other subsurface processes fall to the right of the flow line (consequently, as will be seen from Figs. 3 and 4, drainage basin development is fundamentally a subsurface process) and all atmospheric processes to the left. Thus, the other subsurface velocities indicated on the diagram are less than the subsurface flow rate and all surface velocities are greater. Further, the radial pattern of icons diverging from the upper left suggests that subsurface process rates, except that of water flow, tend to decrease with increasing scale, while flow rates associated with the atmosphere tend to increase with scale until limits of planetary size are approached. The decrease in subsurface velocities with increasing scale is associated with the heterogeneous networks over which the processes evolve and the diminution in connectivity with increasing scale. The increase in velocity of atmospheric processes with increasing spatial scale is a sign of their foundation in non-linear dynamics and is more easily seen in the inverse view where momentum transfer to smaller scales through turbulence is terminated by frictional interaction with heterogeneities on the Earth's surface. This produces a decrease in fluid flow rates with decreasing length scales approaching surfaces. Fluvial processes along the surface appear consistent with a velocity that increases slightly with scale and is certainly larger than the subsurface flow rate. According to Skoien et al. (2003), indeed, the appropriate dependence is represented by t=Ax^0.9, which is in accord with the diagram. This result is generally consistent with a non-linear flow phenomenon interacting with surface heterogeneity.

Figure 3 Figure 2.9 from National Research Council (1991) overlain with predictions of the hierarchical growth of two-dimensional optimal flow paths in the subsurface (green and brown) and the associated prediction of the subsurface flow rate (in blue). The subsurface flow rate chosen was 2 m yr⁻¹, equal to the geometric mean of the values used later in Fig. 3 to show the typical range of annual flow rates. The green path uses the tortuosity exponent of standard percolation, 1.13 (applicable to more homogeneous models), and is associated with a Hack's law exponent of 0.565, while the brown path uses the tortuosity exponent from the optimal flow paths in a highly heterogeneous environment, 1.21, which generates a Hack's law exponent of 0.605. The values 0.565 and 0.605 are shown in Hunt (2016) to give an excellent accounting for actual values of Hack's law relating stream length to drainage basin area (between 0.57 and 0.60) as reported in Maritan et al. (1996).

Figure 4 After Hunt (2017b). “Optimal paths” represents a prediction from scaling relationship (Eq. 1) with $s=\mathrm{0.83}=\mathrm{1}/D_{\text{opt}}$ (2D) and a flow rate that ranges from 0.25 to 25 myr⁻¹. “Solute transport” substitutes the backbone exponent, D_b= for D_opt, but it is otherwise unchanged. The flow rates relevant to vegetation growth and soil formation should be smaller than the range given in Bloeschl and Sivapalan (1995). Data corresponding to “optimal paths” include over 6000 results for tree heights (or root lateral spreads) as well as 31 river basin changes due to stream capture or sill overtopping and, finally, the centroids of the icons for “Shallow groundwater circulation”, “Drainage network formation”, and “Development of major river basins” from NRC (1991). Data for “Solute Transport” correspond to soil depths from dozens of studies around the world, as well as new sources listed here, particularly commencing in the middle Mesozoic. Those with ages in the millions of years up typically describe depths of deep tropical weathering, or laterite soils. “Space limit”, a maximum accessible scale, is chosen as the linear dimension of a supercontinent, which we approximated as the square root of today's total continental area (ca. 12 500 000 km²), while “Time limit”, presumed to be 150 Myr, is the length of time since the break-up of Pangaea (Cawood and Hawkesworth, 2015). Thus, the limits are assumed imposed by the resetting of tectonic structures known as the Wilson cycle.

An increase in groundwater flow rate would be expected to increase advective flow limited processes in the subsurface as well. Without such an increase, increased carbon sequestration is likely to be limited, except for mitigation of the loss of sequestration due to (soil) ecosystem damage. In order to speed up the water cycle, groundwater flow rates would have to be increased, or else the limiting effect of groundwater flow speeds would need to be bypassed. Where water flow rates are flux-limited, an increase in precipitation will be reflected in groundwater flow, but where the hydraulic conductivity is too small to accommodate increased fluxes, the increased hydrologic fluxes will mostly tend to increase surface run-off. But bypassing groundwater flow paths through higher velocity surface paths (due to land-use changes) will only allow an acceleration of the water cycle unaccompanied by effects on the carbon cycle, which are mostly associated with chemical weathering and plant growth.

2.2 The NRC (1991) Stommel diagram

The publication Opportunities in the Hydrologic Sciences, known as the Blue Book (National Research Council, 1991), was fundamental to the foundation of the hydrologic sciences program at NSF. Figure 3, presenting the modified Fig. 2.9 (on p. 59) of Garrison Sposito from the book Opportunities in the Hydrologic Sciences, illustrates ranges of hydrologic process scales. We should point out that the original Fig. 2.9 was reused as Fig. 2.2 in the 2001 report Basic Research Opportunities in Earth Science, which figured in the organization of the CZO observatories. In the 2001 publication, the axes were distorted somewhat such that equal intervals do not correspond to equal ratios of time. We reproduce the original figure here as Fig. 3 overlain by predictions derived from subsurface flow and vegetation growth (Hunt, 2017b).

In Fig. 3, the length scale associated with “soil depth” icon clearly does not refer to depth but rather to horizontal length scales; otherwise we could have a 1000 km deep soil in about 10 000 years.

Consider next the three icons: “shallow ground water circulation”, “channel network formation”, and “development of major river basins” that form an approximate line. What could that line represent? Hunt (2017a) compared percolation theoretical expression for tortuosity with implications for stream sinuosity of Hack's (1957) law relating stream length, L, to drainage basin area, A,

$$\begin{array}{}\text{(2)}& L=C(A{)}^p.\end{array}$$

When miles are used as units of length, C=1.4 (Rigon et al., 1996), but for kilometers, C is (1.4)(1.6)=2.24. The values of the exponent p are found to be between 0.57 and 0.6 (Maritan et al., 1996, though Rigon et al., 1996, argue for 0.6), which produce stream sinuosity (ratio of stream to basin length) exponents exactly twice as large (between 1.14 and 1.2) because the basin area basin (not river) length, l, relationship is Euclidean, i.e., $A=C^{\prime }{l}^{\mathrm{2}}$ (Montgomery and Dietrich, 1992). The tortuosity exponents of percolation theory (1.13 for random networks, and D_opt=1.21, for strongly heterogeneous networks, Hunt and Sahimi, 2017) thus described rather precisely the range of observed stream sinuosity. Our question here is whether there is a link through these exponents to the fundamental pore network scale. And does this link then extend to a spatiotemporal scaling function (Eq. 1)?

If we use the geometric mean annual flow rate from Fig. 4 (2 m yr⁻¹) and the same fundamental length scale of 30 µm together with either the simple tortuosity exponent, 1.13 (the green line), or the optimal path exponent, 1.21 (the brown line), the three considered icons relating to shallow groundwater circulation and drainage basin development are covered nicely. We wish to investigate this possibility below. First we check whether evidence exists that groundwater flow can be important to drainage development.

Petroff et al. (2013) related river drainage development to subsurface flow patterns by looking at flow convergence in areas with insignificant relief, such as Florida. The verified relationships extended to bifurcation angles. Brocard et al. (2011) and Yanites et al. (2013) have also noted a potential role of subsurface flow in stream capture in that groundwater piracy may accelerate mass wasting and erosion of an interior divide. Laity and Malin (1986) and Baker et al. (1990) demonstrate that groundwater sapping (flow convergence but also seepage-induced chemical erosion) controls the rate of headward migration of drainages over long periods of time, as well as the architecture of amphitheaters. However, tectonic processes also contribute to changes in subsurface potentiometric gradients. Thus, important roles of subsurface flow in drainage basin development have already been recognized, and the present implication merely extends the apparent relevance to other properties as well as larger length and timescales.

Given the possible influence of subsurface processes on drainage basin development, the seeming correspondence of the icons of the Stommel diagram in Fig. 3 that represent river basin development, and the apparent coincidence of plant root development and river path development along optimal paths in a highly heterogeneous environment, we propose now investigating the spatiotemporal scaling of drainage basins in terms of Eq. (1), which was found to predict the range of root lateral extent values in the subsurface over timescales up to 100 kyr. Our results are shown in Fig. 4. Although the results regarding tree root growth and river integration appear to be closely analogous, for now we restrict our argument to the mathematical analogues of scales, x₀ and t₀, and form (the value of the exponent), as a physical interpretation of this analogy may require discussion in greater depth and further investigations.

Before discussing the interpretation of Fig. 4, we discuss the variability in the theoretical predictions. For plant roots, this variability is represented in terms of variable flow rates, while the variability of the fundamental network scale is smaller and is neglected. For drainage basin development, we will assume almost the same variability in the parameters that define the fundamental subsurface variability, namely a typical network scale of 30 µm as well as the variability in flow rates (between 0.25 and 25 m yr⁻¹).

These results for drainage reorganization can be quantified in terms of length scale (full length of captured river or tributary) and time. The latter is measured starting with (typically) a trigger event from tectonics (sometimes a previous stream capture) and ending with the integration of the drainage. Such length and timescales, while not arbitrary, are also not unambiguous, and considerable uncertainty is present. Frequently, neither the date of the initiation of the drainage basin change nor that of its equilibration is known accurately. For dates that were given in geologic terms, such as Early Pleistocene, we used the standard mean of the stated interval. Dates of tectonic triggers are more broadly defined, while discrete steps of drainage integration (or disintegration) due to steam captures are sometimes more narrowly defined. Sometimes the authors did not give river lengths but only drainage basin areas. In such cases, Eq. (2), Hack's law, could be applied to estimate a river length.

2.2.1 Multiple events: stream captures

Stream capture occurs by headward erosion of, typically, a second stream at a lower elevation or with greater flow. Such progressions are considered to be bottom-up processes. Here, four cases are discussed in which two stream capture events yield two experimental points. The capture of the upper Amazon by the lower Amazon (total length 6400 km) occurred at about 10 Myr before present ca. at least 56 Myr since tectonics began to impact an established Late Cretaceous flow in nearly the opposite direction (Hoorn et al., 1996; Filgueiredo et al., 2009; Albert et al., 2018). In response to this capture event which increased trunk flow, the Amazon basin began to expand northward (Hoorn et al., 1995) with the result that the Rio Casiquiare, a tributary of the Rio Negro, itself a tributary of the Amazon, is currently capturing the Orinoco River. The Casiquiare (356 km) and upper Orinoco (ca. 480 km) combined are about 840 km long, and the time required for integration is on the order of 10 Myr. Information on related drainage reversals was inadequate for the present purposes. A second such sequence of events started with the capture of the upper Colorado (overall length 2300 km) by the lower Colorado about 11 Myr after the split of the lower Great Basin (17 Myr before present, Young and Spamer, 2001) at about 6 Myr before present. Partly in response to that capture, a 290 km long tributary of the Colorado (the Gunnison) changed course around 1 Myr before present (Aslan et al., 2014) and abandoned Unaweep Canyon. A third such sequence of stream captures in the Appalachian Mountains is discussed in Prince et al. (2011), with a 225 km² basin connected to an existing stream in 1–2 Myr, and a 7000 km² basin expected to be captured in a total of 5 Myr. For these two examples, Hack's law (initially derived for this same physiographic province) had to be used to convert basin area to stream length. Yanites et al. (2013) discuss in their Fig. 2 two steps of the reorganization of the upper Rhine River: a 1.3 Myr duration capture of the Aare/Doubs (ca. 450 km) and then a 1.2 Myr later capture of the Rhine above Waldshut (300 km).

In two cases, because dates of original triggering events are not available sequential captures yield only a single data point. Along the course of the Meuse, Benaichouche et al. (2016) state “the oldest piracy (leaving the Bar river in an oversized valley) corresponds to the capture of the Haute Aisne by a tributary of the Oise River. The time of the piracy was ∼900 kyr ago. Later, at ∼250 kyr (and about 150 km upstream) the Ornain, and Saulx rivers were captured by the Marne River”. The length scale of 150 km thus corresponds to the timescale of 650 kyr. Each capture shortened the original Meuse River drainage. In Guatemala Brocard et al. (2012) summarize their findings: “these diversions occurred during the last 600 kyr with the following succession of events: (1) capture of the Cahabón River near Santa Barbara at ∼500 kyr; (2) capture of the Cahabón River near Purulhá (ca. 50 km downstream) sometime later (240–450 kyr ago)” (mean of 345 kyr ago). The case for using the particular time and space intervals here is not as clear as with the Meuse (Bernaichouche, 2016), for which the succeeding captures were developed from different streams whose drainages were subparallel and integrated. In the case of the Cahabón River, two unrelated drainages were involved in the captures. Both were considered from the perspective of the disaggregation of a river drainage.

2.2.2 Multiple events: sill overtopping

Sills are high ground between disconnected basins. If a river flowing into a basin fills that basin with a combination of water and sediment, then it may flow over the impediment. The integration of drainages by sill overtopping is a top-down process, distinct from headward erosion. The initiation of the Mojave River drainage system in California is considered to have occurred about 3.5 Myr ago, and it became integrated to a length of 200 km in the last 20 kyr. The Mojave River was dammed at Afton for 160 kyr through pluvial climates before it finally breached the sill and advanced to the Soda Lake about 40 km downstream (Reheis et al., 2012). Less than 10 kyr later, the Mojave River arrived in Dumont lake, 50 km further, and likely reached Death Valley, nearly 150 km further downstream, but Holocene drying interrupted integration of this drainage system (Enzel et al., 2013). If the Pleistocene pluvial period had continued, its full integration to the distance of Death Valley would have been possible in the future.

Per Larson et al. (2020), “A ca. 2.5 Myr [before present] date for the initiation of top-down integration of the Verde River from the upper Verde Valley into what are now downstream basins is consistent with the presence of a 3.3 Myr [ago] volcanic tephra…The basins depicted here were formerly endorheic, but integrated within the last ∼2.2–2.8 Myr. The integration of these basins resulted in the modern through-flowing drainage networks of the (320 km) Salt, (272 km) Verde, and Gila Rivers of central Arizona”. The same time frame was implicitly extended to the Salt River, but not necessarily to the Gila River, of which the Salt and the Verde are both tributaries.

2.2.3 Parallel capture events: durations of events assumed similar

Pastor et al. (2012) examined 14 (parallel but not simultaneous) stream capture events in the region south of the High Atlas in Morocco which occurred over roughly 100 kyr time intervals. The lengths of the drainages involved ranged mostly from 20 to 40 km, from which we extract two data points. Isotopic evidence from the Indus delta indicates the capture of four Punjabi rivers at about 5 Myr ago (Clift and Bluzstajn, 2005). The authors state, “exhumation histories for the western Greater Himalayas shows that these mountains were in existence before 20 Myr ago”, indicating a time period of approximately 15 Myr for their capture. The rivers are between 500 and 1400 km in length and are named Ravi, Sutlej, Chenab, and Jhelum.

2.2.4 Single events

Stokes and Mather (2003) addressed the connections of the east-flowing 100 km long Almanzora River across the Sierra Almagro. Quoting the authors: “Emergence and establishment of continental conditions within the basins to the north and south of the Sierra Almagro during Salmerón Formation time (Late Pliocene–Early Pleistocene [midpoint 2.56 Myr ago]). The connection probably took place at some point between the Plio-Pleistocene and Pleistocene stages of drainage evolution (Early–?Mid-Pleistocene)” (1.43 Myr ago). The time difference is roughly 1.1 Myr. The 4180 km long Niger River drainage began to form between 45 and 40 Myr ago when the inland part of its drainage began to emerge from a shallow sea and connection of the upper and lower reaches was established between 34 and 29 Myr ago (Chardon et al., 2016), yielding 11 Myr. Struth et al. (2020) discuss the reorganization of the 170 km Suárez River basin, including its piracy of additional smaller basins along its east side, over a 405 kyr period from its capture by the Magdalena. However, specific distances for the smaller events were not possible to extract. Goudie (2005) describes the reversal of the eastern 220 km long portion of the proto-Katonga River due to the westward migration of its headwaters in the “swamp divide” over the ca. 400 000 years (Johnson et al., 2000) since the formation of Lake Victoria through downwarping. Fan et al. (2018) demonstrate the reorganization of the Daotang basin within 80 kyr of the capture of Yihe River by the Chaiwen, adding 25 km² to the Yihe River drainage. Hack's law was used to generate a distance from the basin area.

The Yellow River reorganization included a possible combined capture at the edge of the Tibetan Plateau near the Gonghe Basin and top-down organization involving filling of basins from upstream (Craddock et al., 2010; Harkins et al., 2007). Late Miocene to Early Pleistocene sediment accumulated almost continuously in basins until the integration of the Yellow River at 1.8–1.7 Myr ago (Zhang et al., 2014). The reorganization was likely triggered due to initiation of faulting in the range 10–11 Myr ago (Meng et al., 2020). The length of the Yellow River above the Gonghe Basin is about 1300 km. The approximately 100 km long paleo-Daotanghe (which means reversed river in Chinese) was cut off from the Yellow River 500 kyr later by isostatic adjustment to the Yellow River-triggered denudation, with its modern remnant of 40 + km flowing in the opposite direction to a closed basin (Zhang et al., 2014).

2.2.5 Disputed timescales

In their Fig. 6a–c, Wang et al. (2018) demonstrate a 50 Myr duration of the reorganization of the 6400 km long Yangtze across the Jianghan Basin, including flow reversal for the upper half. It should, however, also be mentioned that Su et al. (2019) place the capture of the middle Yangtze from the Red River by the lower Yangtze at 5 Myr ago, 45 million years after its alleged integration by 50 Myr ago (Wang et al., 2018). We used the longer interval.

The integration of the Blue Nile drainage with its length of 1450 km is relatively well constrained from the emplacement of the volcanics between 30 and 29 Myr ago that capped the Ethiopian plateau to captures at 10 Myr ago and at 5 Myr ago that led to 2-fold and 5-fold increases, respectively, in its sediment load (Giachetta and Willett, 2018). (However, Fielding et al., 2018, conducted a detailed geochemical investigation of the Nile delta to greater depths than previously, producing a more nuanced interpretation including much earlier trans-Saharan connections.) We used the longer interval, since the more detailed evidence of the later study does not controvert the specific conclusions of the earlier one.

Suhail et al. (2020) suggest the integration of 500 km of the upper Dadu river in the last 3.8 Myr, although since Yang et al. (2020) suggest that the reorganization occurred within 1.4 Myr, or even 0.6 Myr, the point at 1.4 Myr is also plotted together with the suggested 3.8 Myr.

To confirm the validity of Eq. (1), we performed a regression analysis of lengths predicted from Eq. (1) using v₀=25 m yr⁻¹ versus the observed values, summarized in Table 1. A statistical analysis yields the relationship given by predicted length $=\mathrm{1.081}\times$ integrated length, with an estimated standard error of 8.3 % and an R² value of 0.85. (Note that Table 1 includes 31 data points, but the regression analysis was performed using 29 data points, as two data points – for the Niger and the Amazon – were identified as outliers and excluded from the regression analysis.) Out of 29 analyzed data points, 27 points fall within the 95 % prediction interval, confirming the validity of Eq. (1).

Table 1 Comparison of predicted and actual river lengths.

A search for subsurface flow rate scaling information turned up an excellent suite of nine measured basin residence times for deep groundwater as a function of length scale as measured by ⁴He fluxes from South America (Aggarwal et al., 2014). That study does not support a scale-independent flow rate; rather it supports a decreasing flow rate with increasing length scale. In fact, the observations are in almost exact agreement with Eq. (1) for optimal path flow as the measured exponent is 0.824, less than 0.4 % different from $\mathrm{1}/D_{\text{opt}}$, while the numerical constant implies a flow rate of 35 m yr⁻¹ at a length scale of 30 µm, only ca. 50 % larger than the assumption in Hunt (2017b) and reapplied here (25 m yr⁻¹). Furthermore, the R² value is 0.84. Adding these values to the analysis of the suite of river lengths investigated above leaves the discrepancy unaltered at 8 % and changes R² only minimally.

3 Uncertainties

In addition to uncertainties in measurements and model parameters common to all the processes discussed, drainage basin evolution adds the difficulty of dating and measuring a complex process operating in three spatial dimensions, for which much of the immediate evidence is destroyed over time. A motif that is repeated in the literature discussions referenced above is a kind of dichotomy between the perspectives of (1) a new establishment of a drainage and (2) the evolution of an existing drainage, which may in many ways resemble its present form. Such evolution will typically contain events of tectonic and climatic origins, both leaving some imprint on the hydrology, whose evolution may, at times, be discrete and datable. Thus, what from a static perspective looks like a sudden drainage basin connection appears, in a more integrated perspective, often to be geographic shifts of an existing drainage through smaller events. Geomorphologists, hydrologists, sedimentologists, geochemists, and geodynamicists have differing perspectives, but in fact drainages may have existed a long time already, affecting both connection times and connection length scales, and increasing the number and complexity of the shifts. A wide range of additional uncertainties in interpretation exists, which will not be discussed in detail here, but it is clearly to be expected that the dates and spatial scales of fluvial evolution will themselves evolve in the next decades. If inferred connectivity timescales turn out to be much smaller than generally accessed here, that evidence will favor a greater role of top-down drainage evolution (Hilgendorf et al.,  2020), by basin overtopping, than by headward erosion. The former proceeds more importantly over the surface than through subsurface hydrologic connections, which influence headward erosion. If the overtopping interpretation is correct, we might expect that the evolution of drainage basins would ultimately more nearly follow the results for fluvial process icons on the Stommel diagram of Fig. 2. Although the cases quoted here from the arid southwest fit well with the timescales for headward erosion from more humid climates, large changes in climate from the Pleistocene glacial episodes to the (short) interglacials may complicate interpretation that subsurface flow rates may be significant even when organization of strings of endorheic basins proceeds by sill overtopping.

Consider the sensitivity of Eq. (1) with s=0.83, to its physical parameters, the network scale, x₀, and the flow rate, v₀. Representing t₀ as $x_{\mathrm{0}}/v_{\mathrm{0}}$ clarifies that the length scale enters as $x_{\mathrm{0}}^{\mathrm{0.17}}$. Three orders of magnitude change in x₀ changes x only by a factor 3. In contrast, x is proportional to $v_{\mathrm{0}}^{\mathrm{0.83}}$, a much higher sensitivity. Division of integrated drainage lengths by known values of $v_{\mathrm{0}}^{\mathrm{0.83}}$ would allow universal graphical representation, aside from the slowly varying factor $x_{\mathrm{0}}^{\mathrm{0.17}}$. Using a typical global value for v₀ in Eq. (1), while necessary at this time, neglects significant variability in flow rate, a likely source for much of the statistical variability in the regression.

4 Limits

Consider Fig. 4 once more. The maximum timescale found for soil depths is 150 Myr, and river basin organization (100 Myr) corresponds to mid to late Mesozoic, if continuous development to the present is appropriate. The predicted maximum drainage basin length at 150 Myr is ca. 10 000 km, not very different from the size of a supercontinent and not so different from the timescales proposed for reorganization of the ca. 6400 km long the Amazon and Yangtze (50–100 Myr). The middle of the predicted range of soil depths is a couple hundred meters, also corresponding to observation of the vertical extent of the critical zone. Notably, 150 Myr before present corresponds also to the breakup of Pangaea (Cawood and Hawkesworth, 2015). As stated by Gupta (2007) “There is thus a link between the evolution of large river systems with the Wilson Cycle – the creation and destruction of supercontinents, such as Rodinia and Pangaea”. Given the geologic reworking of Earth's surface (and reorganization of mountain ranges as well as oceanic boundaries) on a Wilson cycle timescale, we expect a significant drop-off in the availability of relevant and uninterrupted hydrologic and soil depth data at or near the timescale since the breakup of Pangaea. This discussion justifies the “limits” in time and space put on Fig. 4.

We can use the predictions that underlie Fig. 4 to generate an exact (though not necessarily accurate) prediction of the architecture of basins at 150 Myr, namely a maximum weathered profile depth and a maximum linear extent (Table 2). We apply the maximum and minimum flow rates from Fig. 4 as well as the geometric mean. Use of the same flow rates for soil formation (infiltration fluxes) and vegetation growth (transpiration fluxes) is reasonable, even though one can expect the latter, on average, to be about a factor of 2 larger (Hunt et al., 2020b). Some uncertainties may arise, however, from using the same range of flow rates for all three phenomena simultaneously. In spite of these complications, which may lead to a modest overestimate of soil depth, we simply use the same velocities and the same fundamental length scale that led to the predictions of Fig. 4. Clearly, these values are on the right order of magnitude. The estimations may help constrain sediment and carbon fluxes over geologic timescales.

Table 2 Predicted spatial scales at 150 Myr for various flow rates.

Consider the dependence of maximum drainage basin length on flow rates. The longest interconnected drainage under arid to semi-arid conditions (P on the order of 10 cm yr⁻¹) is predicted to be only 248 km. Although this may be a bit extreme, it is also clear that in truly arid regions of under 10 cm yr⁻¹ (central Australia, Saudi Arabia, and the Sahara) development of connected long drainages is uncommon. As pointed out by Maxwell et al. (2016), subsurface flow rates derived from residence times are proportional to the difference of precipitation and evapotranspiration. Thus, while wet climates can generate a drainage basin of continental size (over 10 000 km), in dry climates, lack of connectivity of drainage basins likely prevents the organization of such large basins. Further, drying climates will disconnect and inactivate large drainages, such as the drying of the Sahara since the end of the Miocene (Goudie, 2005) or even, in the Holocene, the incipient Pleistocene Mojave River connections (Reheis et al., 2012; Enzel et al., 2013).

As for specific depths of weathered profiles dating back to Mesozoic times, values of about 130 m have been reported in Germany of near 150 Myr (Felix-Henningsen, 2018), of 150 m near Belgium (Felix-Henningsen, 1994), of 120 m in Australia (though not continuously forming over the entire interval; Gardner, 1957), over 200 m in southwestern England (Migon and Bergstrom, 2001), and in excess of 100 m in South America (Giovanni et al., 2017; Ruckmick, 1963), with the deepest at 240 m (Ruckmick, 1963). Tardy and Roquin (1992) state: “Over most of the areas of the Brazilian and the African shields, a very thick lateritic mantle has been continuously formed over more than 100 Myr. Laterites are widespread in Australia, India, Burma, Brazil and in intertropical Africa. In these regions, bedrock is generally weathered to depths of over 10 to 150 m”. Other references to very old weathered profiles are given in Yu and Hunt (2017a) and Hunt et al. (2021).

Note that in Fig. 4, erosion is not considered, though it has been addressed in Yu and Hunt (2017a, b) and Egli et al. (2018) in greater detail. Its full treatment is more complicated than given here. However, it is known that erosion rates as small as 1 m Myr⁻¹ occur only in arid continental interiors (Bierman and Nichols, 2004), such as Australia, where the precipitation can be as low as 0.04 m yr⁻¹. In wetter climates, even a denudation rate of 5 m Myr⁻¹ would add up to 750 m over 150 Myr timescales, which means that our estimates of weathering profile depths are probably not particularly high.

Figure 5 A demonstration that tectonic (constant velocity) and predicted optimal path scaling, when using space and timescales defined at the pore scale and for pore-scale flow rates, converge at timescales larger than 1 Myr and length scales exceeding 100 km. In this, consider that 25 m yr⁻¹ is less than a maximum flow rate (35 m yr⁻¹ was also found here from data of Aggarwal et al., 2014) and that slower flow rates may also be relevant. For comparison see “Optimal (2.5 m yr⁻¹)”, which shows predictions of Eq. (1) with a flow rate of 2.5 m yr⁻¹ instead of 25 m yr⁻¹. If, as appears possible from the figure, drainage basin scaling results are found in the wedge of time and length scales between optimal path scaling and tectonic scaling, the range of drainage basin sizes becomes restricted more sharply due to limited water supply at shorter timescales than the values defined by the intersection of optimal paths with the Wilson cycle timescale. In particular, the timescale at this intersection for a flow rate of 1 m yr⁻¹ and a tectonic rate of 0.035 m yr⁻¹ is about 300 kyr and the length scale is about 13 km.

Interaction with spatiotemporal scaling associated with tectonic drivers may introduce an important cross-over in the variability of drainage basins with scale. According to Roberts (2019), “At large scales (≳10 km, ≳1 Ma) drainage networks appear to have a synchronized response to uplift and erosional processes. At smaller scales erosion generates complex landforms. Most power and commonalities reside at long wavelengths and timescales (>100 km, >1 Ma). The presence of commonalities between drainage networks is expected for systems in which large signals (e.g., uplift) are forced through random media (e.g., lithology, biota)”. Such an increase in commonality may reflect an interaction between relief producing (tectonic) and relief shaping (drainage basin organization) processes, as Roberts (2019) says, “a synchronized response”. Consider Fig. 5, where only the river drainage data are retained. Typical horizontal tectonic velocities are ca. 2–3 cm yr⁻¹ (though oceanic plates often reach 6 cm yr⁻¹), a time-independent velocity. The reference to the key length and timescales, 100 km and 1 Ma, generates an upper limit on horizontal tectonic velocities of ca. 10 cm yr⁻¹. In Fig. 5, tectonic rates of 3.5 and 10 cm yr⁻¹ are given, together with the groundwater flow rate and its associated optimal path scaling function from Eq. (1) using the same parameters as in Fig. 4 but including a range of flow rates of 1 order of magnitude. Because the slopes of the scaling relationships are so similar (1 for flow, 0.83 for optimal paths), the relatively small variability in parameters can change the intersection of these lines from 1 kyr to 100 Myr. The maximum tectonic rate and the given 25 m yr⁻¹ flow rate, for example, meet in the upper corner of the diagram.

The data found for drainage basin reorganization commence at about 100 kyr and 10 km. Visual examination of Fig. 5 may suggest that drainage basin development at longer timescales is confined between the tectonic and the optimal path scaling functions. A related implication seems to be that, when the optimal path function falls below the tectonic relationship, no timescale will suffice for equilibration of hydrology to tectonic changes that plot up beyond that timescale. This restriction may be especially important in arid regions.

Use of a groundwater flow rate of 10 m yr⁻¹, just over a factor of 2 smaller than applied, leads to equality of tectonic and river spatial scales of 100 km at exactly 1 Ma. If we had used the fundamental groundwater flow rate extracted from Aggarwal et al. (2014) of 35 m yr⁻¹, virtually all the drainage basin data would lie between the optimal path scaling and the typical tectonic rates. Finally, for intermediate values of these fundamental velocity parameters, there is a large range of timescales starting near 1 Myr and extending to the end of the diagram, in which the time and length scales of Eq. (1) and the optimal path scaling match up with those from tectonics. This approximate equality may help explain Roberts' (2019) “synchronized response” from above.

In his opening sentences, Roberts (2019) states: “Evolution of Earth's surface is a result of geological and geomorphological processes operating at a range of spatial and temporal scales. Therefore, a framework that unifies short wavelength (10 km) process-orientated landscape evolution models with longer wavelength observations and phenomenological approaches is attractive”. One should then conclude that a theoretical treatment of optimal flow paths unifying observed space and timescales at the pore scale with observed vegetation growth at scales up to 100 kyr and 10 km, and with continental wavelength observations at up to the Wilson cycle for rivers, is interesting as well. That the solute transport scaling unifies the pore scale and the weathered depths over the same range in timescales is probably of additional interest.

5 Summary

We have introduced predictive scaling equations in subsurface hydrology that incorporate fundamental and universal spatiotemporal scaling relationships from percolation theory populated by parameters that represent the fundamental nature of the porous medium as a network and flow rates as given by hydrologic variables. These individual relationships were confirmed separately by comparison with a great deal of data that examined their individual predictions for soil depths and vegetation growth. Here, the implications of these (non-linear) scaling equations to the general organization of hydrologic processes with their roots in subsurface flow are considered and some processes not necessarily recognized as being limited by subsurface flow also identified, such as river drainage organization. In the context of a wider understanding of hydrology, it is important to keep in mind that non-linear spatiotemporal scaling does not necessarily imply non-linear dynamics.

Non-linear dynamics and scaling in the atmosphere and in fluvial processes inevitably traces back to the applicability of the Navier–Stokes (NS) equation. Turbulence in fluvial and atmospheric systems is the rule but laminar flow the exception. In contrast, Reynolds' numbers for groundwater flow are typically on the order of 10⁻⁴ (Hunt and Manzoni, 2015); thus, the momentum transfer in the NS equation is ordinarily negligible. Turbulent flow bounded by a wall (e.g., a river, or the atmosphere, by the Earth's surface) conforms to the “law of the wall”, an increasing mean velocity with distance from the surface (Cole, 1956). Turbulent flow also exhibits an energy cascade associated with a well-characterized reduction in power with decreasing length scales (Kolmogorov 1941), although the introduction of heterogeneity complicates theory (Li and Meneveau, 2005; Vindel et al., 2008). In river (open channel) flow on the Earth's surface, for a given potential gradient, the depth-averaged velocity increases with overall depth approximately according to Manning's equation (Nezu and Nakagawa, 1993). Thus, dominant effects of static disorder (heterogeneity) in the subsurface and dynamic disorder in the atmosphere (non-linearity), although both describable using non-integral power laws (or logarithmic dependences), lead to fundamentally distinct non-linear scaling; in the former, velocities tend to be non-increasing functions of length scale and in the latter, non-decreasing. Further, in the former, the role of the reproducibility and predictability of experiments is emphasized; in the latter the irreproducibility and sensitivity to initial conditions is paramount. This is not to say that surface heterogeneity is not expressed in the atmosphere, nor that chaos is absent in the subsurface. In fact, a kind of slow chaotic process, that of plate tectonics, operates at more or less constant speeds but produces events, such as rift formation and plate collisions, that appear as sequences of discrete events in the geologic record. Each major event, such as the aggregation of supercontinents, or their breakup, generates a cascade of related events, which lead to drainage reorganizations.

The most important interaction between the subsurface and the atmosphere or rivers is the surface roughness, which has contributions at scales from the particle to the continental. Although individual particle length scales have the dominant impact on subsurface flows, continental scales have the greatest impacts on atmospheric fluxes (e.g., Stendel et al., 2021), while rivers can interact with this roughness at all length scales. We have neglected complications due to hyporheic exchange (Bencala, 2006) of water fluxes between the surface (rivers) and the subsurface within a drainage network, although they will presumably lead to some mixing of surface and subsurface characteristics at the interface.

An important question is, therefore, over how wide a range of length scales can our predicted scaling relationships hold. The depths of soils (or deep tropical weathering) are available for timescales up to about 150 Myr, with the deepest exceeding 200 m. Widespread depths exceeding 100 m are cited. The time for river basin organization to reach continental scales, on the order of 10 000 km, also appears to be about 150 Myr. Both of these numbers are in accord with the prediction of the scaling relationships accessed. That the two already existing scaling relationships predict reasonably accurately the largest horizontal and vertical scales of what is essentially the critical zone over a geologic timescale equal to that since the separation of Pangaea suggests anew the fundamental relevance of the Wilson tectonic cycle in constraining our knowledge of hydrogeomorphic evolution. The coincidental near equivalence of the optimal path scaling function and its slow decline of velocity over time with the velocities associated with tectonics at timescales from approximately 1 to 100 Myr seems likely to form the basis for the commonalities of drainage basin evolution displayed on these timescales, as well as the greater variability at shorter timescales, particularly on smaller regions of continental land mass. It also suggests a greater relevance to continental-scale processes of flow and transport paths traced all the way back to pore scales than hitherto typically assumed. Subsurface flow paths and rates appear to be critical to drainage basin development.

Our concluding points follow:

• Approaching spatiotemporal scaling diagrams with a network theory of flow and transport allows spatial and temporal scales of seemingly different processes to be unified in terms of concepts, theories, and flow rates.

• The scaling relationships can hold over timescales from seconds to hundreds of millions of years.

• Hydrologic observations appear to be constrained significantly by the Wilson tectonic cycle of supercontinent origin and breakup (most recently, 150 Myr).

• Figure 4 shows that the subsurface water flow rate is a critical control of hydrologic processes, limiting both managed (crops) and natural plant growth, as well as soil formation, and, in this preliminary investigation, appearing to limit river drainage development.

• Figure 5 shows the near coincidence of tectonic and optimal path length scales over a wide range of timescales beyond 1 Ma, perhaps providing a mechanism for more universal drainage basin characteristics in this range of timescales.

• Relevant processes for which spatiotemporal scaling reveals characteristic velocities slower than the subsurface flow velocity appear also to have velocities that diminish with length scales. As such, we suggest that their predominant influence is from heterogeneity.

• Processes for which spatiotemporal scaling reveals characteristic velocities faster than the subsurface flow velocity appear to be associated with velocities that increase with length scales. We suggest that these processes are primarily influenced by non-linear dynamics. Such velocities can reach a maximum on spatial scales of continental separations.

• Clearly, surface flow is influenced by a wide range of surface heterogeneities as well as non-linear dynamics., Fluvial processes (Fig. 1) appear to exhibit scaling with a characteristic velocity slightly higher than subsurface flow rates and with a velocity that appears to increase slightly with increasing scale. We suggest that this result implies a somewhat greater relevance of non-linear flow equations than heterogeneity to the dynamics of surface flow, in spite of the apparent dominance of heterogeneity in the subsurface to the evolution of the shape and position of the flow paths and the architecture of drainage basins. Such a result is also in accord with suggestions that river drainage reorganization by overtopping of barriers to flow should occur more rapidly than by headward erosion; however, the apparent success in predicting timescales for drainage reorganization using the subsurface optimal flow paths scaling relationship suggests an important role of subsurface flow in some drainage reorganization events, even when the apparent hydrologic controls relate to surface features.

Finally, we hazard a prediction. Willett et al. (2014) use drainage basin disequilibrium calculations to analyze the capture of the Apalachicola River by the Savannah. The portion of the paleo-Apalachicola between the present divide and the Savannah then reversed course. On their map, this distance is roughly 14 km. Use of Hack's law converts this distance to a stream length of ca. 38 km. Using 38 km in Eq. (1) yields a time of 126 kyr for the capture. In future research, we wish to explore further the connection between tree root growth and river drainage organization.