Emergent simplicity despite local complexity in 1 eroding fluvial landscapes 2

Movies S1 and S2 (showing examples of the time dependent behavior of the threshold model), and a simple mathematical explanation for how models of physical erosion can be simplified to few parameters.


INTRODUCTION 19
Landscape evolution is a response to physical, chemical and biologic processes operating 20 across a broad range of scales (e.g. Gasparini et al. 2006, Anderson & Anderson 2010, 21 Egholm et al. 2013). Conversely, it plays an important role in determining geologic, chemical, 22 climatic and biologic evolution, and in distributing natural resources (e.g. Howard et al. 1994, 23 Scheingross et al. 2019, Fernandes et al. 2019. I focus on physical erosion of longitudinal 24 river profiles carved into bedrock, which can set the pace of landscape evolution in low to mid-latitudes (e.g. Young & McDougall 1993, Rosenbloom & Anderson 1994, Howard et al. 26 1994, Sklar & Dietrich 1998, Whipple & Tucker 1999, 2002. 27 28 Many natural systems, including rivers are characterized by scale dependent complexity 29 (e.g. Roberts et al. 2019). At small scales, e.g. < O(10) km, < O(10 3 ) years, erosion rates 30 (and river profile evolution) can be highly variable and dynamical physics-based models 31 have been developed to understand observed complexity (e.g. Lamb et al. 2008). At larger 32 scales river profile evolution is often explored using phenomenological (essentially 33 kinematic) models that capture emergent simplicity (e.g. stream power model; Rosenbloom 34 & Anderson, 1994). Preliminary work examining the spectral content of river profiles 35 suggests that their geometries are scale dependent (e.g. Roberts et al. 2019;Wapenhans 36 et al., 2021). These observations suggest that different physical processes control river 37 profile evolution at the diverse scales of interest (up to O(10 3 ) km). 38

39
Here, I seek an understanding of how scale dependent river geometries emerge and, 40 relatedly, how physics-based erosional models and insights from phenomenological 41 approaches might be formally reconciled. I test three hypotheses. First, most erosional 42 processes can be described using simple threshold models. Secondly, erosional thresholds 43 are responsible for generating complexity at small scales and emergent simplicity. Finally, 44 simple threshold models can replicate predictions of phenomenological (e.g. stream power) 45 models. The approach taken here is akin to particle-based landscape evolution modeling 46 (e.g. Lamb et al. 2008, Tucker & Bradley 2010. 47 yield low residual misfits to observed profiles, despite not explicitly considering processes 51 that determine erosion rates at some scales (e.g. hydrodynamics, substrate changes; Sklar 52 & Dietrich 1998, Tinkler & Wohl 1998, Stock & Montgomery 1999, Lamb & Dietrich 2009, 53 Roberts & White 2010, Salles 2016, Fernandes et al. 2019, Glade et al. 2019, Scheingross 54 et al. 2019. One example is the stream power model, which can be expressed as 55 In this scheme, which is usually presented without accompanying noise (η), erosion rate is 59 set by the velocity of kinematic erosional 'waves' that propagate upstream. The rate at which 60 slopes, ∂z/∂x, propagate is set by the erosional prefactor v, upstream drainage area, A(x), 61 and exponent m as a proxy for discharge. If n ≠ 1 propagating slopes can induce shocks 62 (e.g. steeper slopes can overtake gentler slopes; see Pritchard et al. 2009). Elevation, z, is 63 typically added to profiles by assuming the form of uplift rate, U, which can vary as a function Calculated incision is shown in Figure 1c. The gray curves in Figure 1a and 1c show results 69 for additional monotonic noise (η > 0), which generates a few meters of relief. These results 70 indicate that long wavelength structure can emerge through local complexity in simple 71 phenomenological models of landscape evolution (Roberts et al. 2019). The theory of 72 stochastic partial differential equations gives a basis for understanding why some natural 73 phenomena generate relatively simple structures at large scales despite considerable local 74 complexity (Kardar et al. 1986, Hairer 2014. Pioneering geomorphic work has also related erosional processes operating at small scales to larger-scale landscape evolution (Smith & Bretherton 1972, Birnir et al. 2001, see also Dodds et al. 2000, Rinaldo et al. 2014 where L is the width of the column, and F g , F b , F d and F τ are forces related to mass of the 104 rock column, buoyancy, drag and shear. H is the height of the rock column and h1 is the 105 length over which drag acts upon the exposed part of the column (see Supplementary  106 Material). Surface and body forces depend on densities of water and rock, water discharge, 107 gravitational acceleration and the geometric properties of the landscape (e.g. slopes). Most 108 are measurable at the present-day and reasonable bounds can be estimated a short way 109 back in time (e.g. from gauge station data). This simple scheme can be simplified further if 110 we assume that blocks will topple if their height exceeds a critical value, c, which could be 111 expressed as a function of, for example, F g , F b , F d , F τ (Supplementary Material). In this 112 scheme, which is cast in terms of position along a river, x, and time, t, elevation, z, becomes 113 where ∆ = # % − # %*1 . This model requires a starting condition, i.e. z(x) at t = 0. In the first 117 test, a random uniform distribution of elevations is generated and then ordered by height 118 such that the starting condition is monotonic (Figure 1d The threshold model makes other predictions that might be fruitful avenues for further work. 172 They include autogenic growth of waterfalls, waterfall spacing and height, and fluvial 173 terraces at the crest of waterfalls that are curved as shocks propagate but flat along their base. An initial assessment suggests that these predictions are consistent with some geomorphic observations (e.g. Stucky de Quay et al. 2019, Scheingross et al. 2020; Figure  176 1g). 177

CONCLUSIONS 179
The simple threshold models examined predict slopes that can propagate at different 180 velocities at short length and time scales such that steeper slopes can overtake (and 181 cannibalize) gentler slopes. At longer time scales the model retains local complexity (e.g. 182 changes in slope at a range of length scales). In the examples examined, these local 183 changes in relief are nested within longer wavelength slopes that propagate upstream in a 184 predictable way (e.g. Figure 1g-i). This emergent simplicity gives support for the use of 185 phenomenological models (e.g. stream power) to predict river profile evolution at large 186 scales (e.g. Figure 2b, c, f). The models explored in this paper are somewhat arbitrary and  Supplementary Material for 'Emergent simplicity despite local complexity in eroding fluvial landscapes' Gareth Roberts, gareth.roberts@imperial.ac.uk Department of Earth Science and Engineering, Imperial College London, UK Summary This Supplementary Information contains two movies (details follow) and a simple mathematical explanation for how models of physical erosion can be simplified to very few parameters. The simple (few parameter) model is amenable to a straightforward, computationally inexpensive, exploration of parameter space at much larger scales. For example, Figure 2 shows the results of running the model where the evolution of 10 5 blocks is predicted for 10 5 time steps, which takes 50 s using a 2.6GHz Intel Core i7 processor. Finally, results showing the effect of changing the critical threshold value, c, are given in Figure 3 of this document. Results are described in the main manuscript and the movies help to show the time dependent behaviour.

Movies
Movie 1 shows the time dependent evolution of solutions to Equation (3)

Simplifying a physical model of block toppling
The following describes how physical models of erosion along rivers can be described as a consequence of thresholds. The resultant simple models have very few parameters. In the main manuscript a simple (few parameter) model is explored for insights into the evolution of fluvial landscapes from very small (meter) to large (tens to hundreds of kilometres) scales.
Physical erosion is a consequence of body or surface forces (F ) being sufficiently large that erosional thresholds, c, are exceeded. More formally, in discrete notation, at any position along a river, x, elevation will change as function of time, t, such that where ∆z is change in elevation, which can be set by, for example, the size of the rock mass (e.g. pebble, basalt column, fractured schist) being moved between time t and t + 1. This simple description could be expanded to incorporate, for example, shear stresses or drag and critical thresholds for sliding, saltation, toppling or fracturing. The simple model appears to be a universal description of physical erosion along rivers. This supplementary document shows one way in which a simple physical model of blocks toppling (e.g. Lamb & Dietrich, 2008;Stucky de Quay et al., 2019), which appears to be a reasonably description of fluvial erosion in regions of exposed bed rock, can be reduced to a simple model in which erosion occurs if rock column height exceeds a critical value for toppling (i.e. ∆z > c). It is straightforward (and computationally efficient) to expand this model so that the consequences of local physical erosion for fluvial erosion at much larger scales can be explored. Simplification of other well known erosional models (wear; transport-limited erosion) are also examined.
In the simple scheme explored here, the propensity of columns of rock to topple is estimated as a function of drag, shear stress, rock mass and buoyancy. The force generated by drag on the (unit width) column of rock can be expressed as where ρ w is density of water, C d is the dimensionless drag coefficient, u is water velocity, h 1 is height of the column exposed to flowing water. For reasonable values of parameters (see Table  1) in Equation (1), F d is O(10 3 − 10 6 ) N for a column of unit width. The force generated by shear at the top of the unit width column can be expressed as where g is gravitational acceleration, h 2 is depth of the flowing water, dz/dx is channel bed slope, and L is width of the column. F τ is expected to be O(10 − 10 3 ) N for slopes between O(10 −3 − 10 −2 ). The buoyancy force generated as a result of water displaced by the column of unit width rock can be expressed as where h 3 is depth of the water at the base of the column. F b is expected to be up to O(10 5 ) N. The force exerted by the column of unit width rock is where ρ r is density of the rock column. F g is expected to be up to O(10 5 ) N.
Calculating moments (see Figure 1) generated by application of these forces indicates that the column of rock will topple if Substituting Equation (4) into (6) and rearranging to make column height the subject yields If 2F τ + 2F d ≥ L 2 ρ r g, the column will topple if, If 2F τ + 2F d < L 2 ρ r g, the column will topple if,   (b) Power spectrum (from Fast Fourier Transform) of elevation (red circles) used to generate the random noise in the starting condition and relief (black circles). Note elevation spatial series has a white noise spectrum (solid red line), consistent with short wavelength ( 100 km) spectra of some real rivers (Roberts et al. 2019;Wapenhans et al., 2021). Black solid line = power ∝ k 2 , where k is wavenumber. (c) Histogram showing distributions of elevations (red) and relief in the starting condition (black). (d) 100-km-long river profile, containing 10 5 (1 m wide) blocks, evolving for 10 5 time steps. Thick black line = starting condition, thin lines = predicted profile every 10 4 time steps. Threshold, c = 0.5 m in this example. If block toppling occurs at a rate of 1 /year to 1 /century this model represents 10 5 to 10 7 years of evolution.