Scale‐Dependent Flow Directions of Rivers and the Importance of Subplate Support

Large rivers play crucial roles in determining locations of civilization, biodiversity, and efflux to the oceans. The paths they take across Earth's surface vary with scale. At long‐wavelengths rivers can have simple flow paths. At smaller scales, in meanders for example, their paths change rapidly as a consequence of lithology, biota, and other environmental variables. It is not straightforward to identify the scales at which river planforms are set. We overcome these issues by developing a spectral (wavelet) methodology to map flow‐directions as a function of distance and scale. This methodology allows short‐wavelength features (e.g., meanders) to be filtered from river flow‐paths. With short‐wavelength structure removed, the flow‐directions of rivers in Western USA correlate with long‐wavelength gravity anomalies suggesting control by subplate support. This relationship is replicated by an ensemble of landscape evolution models. These results combined suggest that drainage at large scales, O(103) km, is set by subplate support.

While it may seem trivial that rivers flow down topographic gradients, they can follow circuitous routes at a range of scales. At short (<100 km) length scales flow-directions can be extremely variable, which can result in rivers flowing in the opposite direction to the long wavelength topographic gradient (e.g., Goosenecks, San Juan River, USA). These simple observations indicate that river planforms are scale-dependent.
In many instances visual inspection of drainage planforms provides most of the information we need. For example, the Colorado River, which drains western North America, flows mostly to the west, southwest and south in its upper-, mid-, and lower reaches, respectively ( Figure 1). The planform of rivers on top of topographic swells in other continents also have similarly simple patterns at long, O(10 3 -10 4 ) km, wavelengths (Rudge et al., 2015). At these long wavelengths rivers mostly flow away from crests of topographic swells that are supported by the mantle Faccenna et al., 2019;Roberts et al., 2012). This pattern of emergent simplicity at long wavelengths is manifest in the flow paths of many large rivers draining topographic swells and tectonic topography on Earth, for example, African swells, Colorado Plateau, Mexican Highlands, East Australian Highlands, Himalayas, and elsewhere, for example, Tharsis Rise, Mars (Black et al., 2017). These diverse examples suggest that the subplate mantle plays an important role in setting drainage planforms (e.g., Braun, 2010). This contribution has two primary goals. First, we seek to formalize the observations of scale-dependency by developing a spectral methodology to map planforms and flow directions as a function of scale and position. We explore one way in which this approach can be used to compare drainage patterns to environmental variables at appropriate scales. Second, we explore how planforms can be predicted a priori given forcing by subplate support. We achieve this second goal using an ensemble of landscape evolution models.

Wavelet Analysis
A number of approaches could be used to analyze the variability of flow-directions at different scales. Fractal approaches have been used to investigate scales of self-similarity in river planforms (Rinaldo et al., 1993). Simple filtering (e.g., box-car, Gaussian) could also be applied. A spherical harmonic approach has been used to compare flow directions of rivers to long wavelength topography (Black et al., 2017). Spec-LIPP AND ROBERTS 10.1029/2020GL091107 2 of 12 Figure 1. North American drainage overlain on topography and long wavelength free-air gravity. (a) Selected North American drainage networks extracted from ASTER DEM overlain on ETOPO1 DEM. Major rivers in this study are shown by thick black annotated curves; Miss = Mississippi. (b) Drainage overlain on long wavelength (>800 km) GRACE free-air gravity anomalies; contour interval = 10 mGal (Tapley et al., 2005). If admittance Z ∼ 25 mGal/km, calculated dynamic support of western North America is up to ∼1.5 km. Note broadly radial drainage patterns that flow away from crests of positive dynamic topography (red contours). tral approaches have a number of benefits including mapping of spectral power (e.g., red, pink, or white noise), which can be diagnostic of scaling regimes and specific processes. We avoid using standard Fourier transforms because flow directions are not necessarily periodic (i.e., they are nonstationary signals). Instead we make use of continuous wavelet transforms (Daubechies, 1990;Farge, 1992;Kumar & Foufoula-Georgiou, 1997;Torrence & Compo, 1998). There is a precedent for transforming directional time series into the frequency and frequency-distance domains in the atmospheric and oceanic sciences (Donelan et al., 1985(Donelan et al., , 1996. Wavelet analysis has previously been applied to the curvature of river meanders (e.g., Gutierrez & Abad, 2014;van Gerven & Hoitink, 2010) but not to the flow-direction.
Recent wavelet spectral analyses of longitudinal river profiles, that is, elevation as a function of distance, z(x), has shown that the shape of large African rivers is mostly determined at wavelengths >100 km where their power spectra, ϕ(k), can be characterized as red noise, that is, ϕ ∝ k −2 , where k is wavenumber . At shorter wavelengths power appears to have a pink noise spectrum, ϕ ∝ k −1 . These observations give a basis for modeling longitudinal river profiles as systems that possess self-similar scaling and deterministic behavior at long wavelengths that emerges despite local complexity. It gives a basis for understanding why at large length scales, O(10 2 -10 3 ) km, river profiles on top of dynamically supported topography (e.g., Bié dome, West Africa) have common shapes (Roberts, 2019). In this study, we develop wavelet spectral techniques to map flow directions of continental-scale drainage patterns as a function of distance and wavelength.

Data and Methods
This section contains, first, a description of data used to extract drainage patterns and second, methodologies to perform wavelet transformations of a series of directions. Software to perform directional wavelet analysis is provided (github.com/alexlipp/directional-wavelets).
Transforming a series of directions (in the form of an azimuth between 0° and 360°) has several simple steps. First, the drainage data set is extracted from the ASTER GDEM, which has a horizontal resolution of ∼30 m, using Esri's D8 (steepest descent) flow routing algorithms. Second, latitudes and longitudes are resampled along flow paths (e.g., rivers) so that they have equidistant spacing, which makes them straightforward to transform into the spectral domain. In the examples used in this paper δx = 2 km. Third, distances and azimuths are calculated along the path. Local (point-to-point) azimuths, θ(x), are extracted using the gmt mapproject algorithm (Wessel et al., 2013). Note that input is expected to be positions along a river with longitudes and latitudes in decimal degrees and resolution of up to a few tens of meters.
Applying wavelet transforms to a series of azimuths is challenging because the functions are discontinuous-at least one pole contains a discontinuity, for example, sin(360°) = sin(0°). To avoid this issue, we represent the azimuthal series in complex form, a(x). Azimuths can be considered as complex numbers of unit magnitude and variable phase, θ. Making use of Euler's formula any azimuth, θ, can be represented as exp(iθ) with real part cos(θ) and imaginary part sin(θ), which correspond to northings and eastings, respectively. The complex series to be wavelet transformed is (1) where θ varies between 0° and 360°. The azimuth series was normalized to zero mean, for example, c a x a x a ( ) ( ) , prior to transformation. The resultant series of complex numbers was transformed to generate where ψ is the mother wavelet function as a function of scale, s. This function is the wavelet which is multiplied with the complex azimuthal series at different scales to generate the transformed series. ψ is an oscillation of finite duration with maximum unit magnitude and zero mean.
In the examples shown in this paper the mother wavelet is a sixth-order derivative of a Gaussian (DOG), with δ j = 0.1, although we recognize other mother wavelets could also be used (see Torrence & Compo, 1998, m = 6 . Scales are converted to Fourier periods (Torrence & Compo, 1998).
Real-valued azimuths (in degrees) as a function of distance and wavenumber can then be calculated as θ(x, k) = ζ180/π mod 360, where mod is the modulus operator, and ζ is the argument of the transformed (complex) series (Equation 2). ζ is computed as Note that the mean of the complex series, a, is added to the reconstructed complex series in this step. The inverse wavelet transformation is simply the sum of the signal in distance-wavenumber space over scales, for the DOG mother wavelet used in this study (Torrence & Compo, 1998). Note that the subscript x denotes a transformed series. The denominator factor (here 1.7379) depends on the mother wavelet used in the transformation. Real-valued azimuths (in degrees) can be generated from a x Filtered azimuth series can now be generated by solving Equations 4 and 5 between the scales of interest. Filtered river planforms can be estimated from these azimuths by forward geodetic transformation, which returns longitudes and latitudes given a starting position (e.g., the head of the river), azimuths and distances. In this case, distances are scaled so that the final calculated position coincides with the actual river mouth. The WGS84 datum was used to perform the transformation. While we consider only river paths in this study, it is straightforward to generalize this approach to other sequential paths or directional data (e.g., a time series of flow velocities and directions).
An alternative more intuitive methodology is to transform eastings and northings generated from the azimuthal series as presented in the supporting information. As expected, this approach gives the same results as transforming the complex form of the signal ( Figure S1; supporting information).
There are two main sources of uncertainty in the wavelet transformation described above. First, there is uncertainty in the position of mapped river planforms. Our approach is limited to scales greater than the spatial resolution of digital elevation data (tens of meters). The fidelity of mapped rivers was assessed by comparison with independent satellite imagery. At the scales of interest (i.e., >2 km) planforms are accurately reproduced away from flat topography and standing water (e.g., lakes). There is also an uncertainty, δθ associated with measuring azimuths from discrete digital elevation data, which is inversely proportional to distance between cells, L, such that sin( ) for simple east-west Euclidean flow paths, which yields an uncertainty of δθ ∼ 0.9° for ASTER data (δx ≈ 30 m) if L = 2 km. If L = 100 km, δθ ∼ 0.02°.
Second, spectral leakage can generate uncertainties in calculated azimuths. A guide to the fidelity of the wavelet transform is the accuracy of reconstituted series (i.e., generated by inverse transformation), which, for the examples in this paper, match the θ(x) series within a few percent in terms of error of the mean (see Figure S1a). Our approach is spatially limited to the resolution of the DEM we use (∼30 m in this study). Additionally, we assume that channels can be represented as a single one-dimensional channel. For our DEM resolution, this is a reasonable assumption excepting very large channels. Finally, we note that while LIPP AND ROBERTS 10.1029/2020GL091107 this work is focused on single-threaded channels, the methodology could be modified to include anastomosing and braided rivers.

Results
We demonstrate our approach by first transforming flow directions of the Colorado River, before analyzing the wider drainage of the Western USA. The Colorado River flows across the Colorado Plateau, through the Grand Canyon, to the Gulf of California (Figure 1). Figure 2 shows the results of transforming the flow direction of the Colorado River into distance-wavenumber space. Figure 2a shows measured azimuths from an evenly resampled (δx = 2 km) data set alongside the filtered series for wavelengths >100 and >1,000 km. Figure 2b shows Colorado River azimuths as a function of distance and wavelength (1/k). In Figure 2c, the azimuths for the full resolution data set are shown as vectors with the observed river superimposed on top. These vectors are spread broadly uniformly between 150° ≤ θ ≤ 300° (see rose diagram in Figure 2c). Note LIPP AND ROBERTS 10.1029/2020GL091107 5 of 12 that the map has been rotated. The filtered >100 and >1,000 km azimuths and their associated rose diagrams are shown in Figures 2d and 2e. These long wavelength flow directions have, as expected, a smaller spread than the full data set. The long wavelength azimuthal series (λ > 1,000 km) quantifies flow paths mapped by eye in the introductory section, that is, flow to west (∼270°), southwest (∼240°), and south (∼190°) in the upper, mid, and lower reaches of the river, respectively (Figure 2f).
The white curve in Figure 2f shows a pseudo-Colorado River path generated using only azimuths at wavelengths >1,000 km and forward geodetic transformation. This calculated flow path reinforces our assertion that most of the long wavelength structure of the Colorado River is set by just two changes in flow direction. Long wavelength flow directions of the other rivers in this study (e.g., Columbia, Mississippi) are overlain on topography and long wavelength free-air gravity anomalies in Figures 3a and 3b.

Influence of Subplate Support
Gravity anomalies, tomographic models, magmatism and isostatic calculations, which include crustal thickness estimates, indicate that western North American topography is principally a consequence of subcrustal support moderated by tectonic and erosional processes (Atwater, 1970;Fernandes et al., 2019;Wernicke, 1985). A guide to the amplitude and wavelength of subplate support is the transfer function (admittance) between long wavelength free-air gravity and topography (McKenzie, 2010). We note that gravity anomalies at spherical harmonic degrees appropriate for this study, O(1,000) km, are particularly sensitive to upper mantle structure (Colli et al., 2016). In western North America the calculated admittance is ∼25 ± 3 mGal/km at wavelengths >1,000 km, which implies that up to ∼1.5 km of western North American topography is supported by the mantle (Stephenson et al., 2014). Figure 1b shows long wavelength freeair gravity from the GRACE data set filtered to extract wavelengths between ∼800 and 2,500 km (Tapley et al., 2005).
By removing the short wavelength (<1,000 km) contributions to flow directions, we can now compare the planform of rivers draining western North America to putative evidence for subplate support at appropriate scales ( Figure 3). The black vectors in Figures 3c-3e show mean flow directions of the first 500 km of major rivers draining the Colorado-Rocky-Mountains plateaux filtered to remove wavelengths <1,000 km. In all cases flow is directed away from the crest of topography centered on Yellowstone and the Rio Grande rift. Figure 3c shows, for the same region, shear-wave velocity at 75 km depth from a recent tomographic model that incorporates data from the USArray experiment (Shen & Ritzwoller, 2016). Lithospheric thicknesses generated by converting a global shear wave tomographic model into temperature using an empirical parameterization are shown in Figure 3d (Ho et al., 2016). Cenozoic magmatism from the NAVDAT database and surficial geology from the GMNA repository are shown in Figure 3e. There is no obvious correlation between most long wavelength azimuths and crustal thicknesses (see Figure S2) or surficial geology (Buehler & Shearer, 2017). An exception is that westernmost Cenozoic magmatism tends to be concentrated atop the crest of the swell above warm asthenosphere (e.g., Figures 3e and S2; Klocking et al., 2018). All of the analyzed drainage networks show long wavelength flow directed away from regions of low shear-wave velocity, positive gravity anomalies, and embayments of thin lithosphere. These observations, combined with admittance between gravity and topography, are strongly suggestive of topography and drainage patterns maintained by mantle convection.
We investigate the role the mantle plays in maintaining flow directions by comparing flows paths to gravity data ( Figure 4). We start by assuming that long wavelength free-air gravity anomalies are a proxy for subplate support. Under the hypothesis of mantle supported topography, the long-wavelength flow-direction of rivers should flow parallel to the line of steepest gravity descent. The direction of steepest descent is calculated from the first derivative of the gravity field, arctan2(g′ x , g′ y ), where g′ x and g′ y are the first derivatives of the gravity field in x and y directions (Figure 4a). Figure 4b shows the position of rivers generated by forward geodetic transformation of the long wavelength components of western North America's rivers. It demonstrates that the long wavelength components of large rivers draining western North America are parallel/subparallel to the direction of g′ along most of their lengths. More than 75% of the offsets between the Colorado River and g′ are <30° in magnitude. Similar results are obtained for the upper 500 km of the other rivers shown in Figure 4b (see Figure 4c). Flow directions of the Colorado River are shown with  (Shen & Ritzwoller, 2016). Black/gray vectors = average/all long wavelength (>1,000 km) azimuths for first 500 km of major rivers shown in panels (a) and (b); gray curves = river planforms. Dashed curve fringes Colorado Plateau. Y/R = Yellowstone/Rio Grande rift. (d) Flow atop CAM2016 lithospheric thickness map, which was generated by converting shear wave velocities into temperature (Ho et al., 2016). (e) Flow atop Cenozoic magmatism (red circles) from NAVDAT inventory and surficial geology from the GMNA data set (www.navdat.org; ngmdb.usgs.gov/gmna/). Legend shows lithologies colored by age; numbers = age in millions of years.
predicted flow paths from long wavelength free-air gravity in Figure 4d. The largest discrepancies in flow directions at long wavelengths are at the uppermost headwaters (distances < 100 km), in the Grand Canyon area (1,000 < distance < 1,500 km) and where the Colorado River meets the Gulf of California. We suggest that these results indicate the importance of subplate support in setting the planform of large rivers draining western North America.
LIPP AND ROBERTS 10.1029/2020GL091107 8 of 12 (a) Red-white-blue color scale = long wavelength free-air gravity anomalies (∼800-2,500 km) from GRACE data set, which are a rough guide to loci of subplate support (Colli et al., 2016;Tapley et al., 2005). Vectors = flow directions (gradients) calculated using gravity data, that is, g′ (see body text). © = Colorado River; gray curves = full resolution drainage. (b) Colored curves = positions of rivers calculated by forward geodetic transformation of long wavelength (>1,000 km) component of azimuths; colors = angular offset between azimuths of long wavelength components of actual river and g′ (see body text). White curves = full resolution drainage. (c) Histogram (binwidth = 30°) showing difference (angular offset) between azimuths calculated from gravity data (g′) and long wavelength flow paths: red = Colorado River; black = other rivers shown in panel (b). (d) Flow direction (azimuths) of Colorado River: gray/dotted/black = full resolution/>1,000 km filter/800-2,500 km filter. Red = predicted flow directions from long wavelength (800-2,500 km) gravity.

Synthetic and Probabilistic Landscapes
A implication of this relationship is that planforms could be predicted using proxies for subplate support, that is, long-wavelength gravity anomalies. However, planforms are proposed to also be sensitive to preexisting structure (i.e., antecedence) or affected by stream capture/piracy (Anderson & Anderson, 2010). For example, barbed drainage on the Colorado Plateau has been interpreted as an indication that the Colorado River once flowed toward the northeast (Lucchitta, 1972).
To explore whether subplate forcing can "overwrite" the effects of preexisting structure we use landscape evolution models. A similar problem has been examined by parameterizing landscape evolution models using uplift predicted from "backwards in time" global convection models Faccenna et al., 2019). Our approach has three steps. First, we generated an initial condition (a surface) that contains up to 100 m of uniformly distributed random noise. This surface can be modified to test the impact of different distributions of antecedent noise. Second, calculated dynamic topography was added to this surface by converting long wavelength free-air gravity anomalies using an admittance of 25 mGal/km (Figure 1b). We note that this surface is similar to cumulative post-Cretaceous uplift mapped using the distribution of marine rocks in western North America, and its crest coincides with loci of Cenozoic basaltic magmatism (Fernandes et al., 2019;Klocking et al., 2018). Simulations of North American landscape evolution parameterized in a similar way have yielded broadly stable Cenozoic drainage planforms (Fernandes et al., 2019).
In the third step, the resultant surface is eroded so that channels form using the well-known stream power formulation of fluvial erosion, which has the form of a nonlinear advective equation where z is elevation, t is time, and A is upstream drainage area (Howard & Kerby, 1983;Rosenbloom & Anderson, 1994;Tinkler & Wohl, 1998;Whipple & Tucker, 1999). Erosional constants v = 0.99 m 0.4 /Myr and m = 0.4 were calibrated for North America using incision measurements in the Grand Canyon (Fernandes et al., 2019). Topographic sinks were filled, then Equation 6 was solved numerically using the "Fastscape" algorithm as implemented in the LandLab package (Barnes et al., 2014;Barnhart et al., 2020;Braun & Willett, 2013;Hobley et al., 2017). Flow-routing was performed with the "D8" algorithm to calculate upstream drainage area. In all of the tests we ran, channels formed during the first time step (by 5 ka).
Figures 5a and 5b show calculated topography and drainage areas at 5 ka extracted from one model. This simple model does not include the coastline, nor transcurrent, extensional or compressional plate motions or complex hydraulic or geomorphic processes (e.g., thresholds). Despite the simple formulation, the calculated positions of major rivers and their mouths are in similar locations to actual large rivers (e.g., Colorado, Columbia, Mississippi tributaries). On smaller scales there is significant short-wavelength complexity, imparted by the initial distribution of noise (Figures 5c and 5d). These results suggest that the long-wavelength planform is controlled by subplate support, but short-wavelength complexity can be generated by other processes. The long wavelength flow directions of large western North American rivers are therefore likely to be no older than the ongoing mantle support. This result is consistent with geomorphic and provenance studies, which indicate that the present day planform of rivers draining western North America has been broadly stable during the last few tens of millions of years (Blum & Pecha, 2014;Fernandes et al., 2019;Galloway et al., 2011).
Finally, to explore this result further, we ran 10,000 simulations where different distributions of uniform antecedent noise (with amplitudes up to 100 m) were inserted into the starting condition. The locations of major channels (defined as a cell with upstream area >10 10 m 2 ) were extracted from each of these runs. The probability of any cell containing a major channel was determined by dividing the total number of times a cell contained a channel by the total number of simulations. The resultant "probabilistic drainage map" is shown in Figure 5e. This map indicates regions where channels occurred more frequently during the 10,000 simulations. As we saw with the single model run, despite the simplicity of these models, positions of predicted rivers with highest probabilities and their mouths overlap with the catchments of large rivers draining western North America (e.g., Colorado, Columbia  scales, O(1-10) km, planforms are moderated by other environmental variables (e.g., lithology, biota, and hydrodynamics).

Conclusions
We have mapped the scale dependence of river networks so that they can be compared to independent observations at appropriate scales (e.g., of subplate support). A continuous wavelet approach was used to transform the complex form of distance-azimuth series into the distance-wavenumber domain. Results indicate that the flow-directions of western North American rivers (e.g., Colorado, Columbia, Mississippi) are principally controlled by the shape of long wavelength O(10 2 -10 3 ) km topography. The correlation between long wavelength gravity anomalies, subplate shear wave velocity anomalies, magmatism and river flow directions indicates the importance of subplate support in maintaining the flow direction of these rivers and the positions of their mouths. Generation of synthetic landscapes further demonstrates that largescale planform structure of rivers is set by the solid earth and moderated by geomorphic and smaller scale processes. We suggest that comparing mapped spectral power of azimuthal series to other environmental variables could give insight into shape and origin of planforms at all scales.

Data Availability Statement
Software to perform the directional wavelet analysis is provided at: github.com/alexlipp/directional-wavelets and archived at the point of submission at doi.org/10.5281/zenodo.4067111. ASTER Topographic data can be accessed at from https://asterweb.jpl.nasa.gov/gdem.asp.