## Abstract

Landscape topography is the expression of the dynamic equilibrium between external forcings (for example, climate and tectonics) and the underlying lithology. The magnitude and spatial arrangement of erosional and depositional fluxes dictate the evolution of landforms during both statistical steady state (SS) and transient state (TS) of major landscape reorganization. For SS landscapes, the common expectation is that any point of the landscape has an equal chance to erode below or above the landscape median erosion rate. We show that this is not the case. Afforded by a unique experimental landscape that provided a detailed space-time recording of erosional fluxes and by defining the so-called E50-area curve, we reveal for the first time that there exists a hierarchical pattern of erosion. Specifically, hillslopes and fluvial channels erode more rapidly than the landscape median erosion rate, whereas intervening parts of the landscape in terms of upstream contributing areas (colluvial regime) erode more slowly. We explain this apparent paradox by documenting the dynamic nature of SS landscapes—landscape locations may transition from being a hillslope to being a valley and then to being a fluvial channel due to ridge migration, channel piracy, and small-scale landscape dynamics through time. Under TS conditions caused by increased precipitation, we show that the E50-area curve drastically changes shape during landscape reorganization. Scale-dependent erosional patterns, as observed in this study, suggest benchmarks in evaluating numerical models and interpreting the variability of sampled erosional rates in field landscapes.

## INTRODUCTION

Landscape topography is sculpted via material fluxes that are controlled by the interplay of different external forcings, such as climate and tectonics, with the underlying lithology (*1*–*6*). Landscapes evolving under constant external forcings tend to achieve steady-state (SS) configurations, where the material flux provided by rock uplift relative to base level is balanced by erosion. These landscapes can be subdivided into different geomorphic process regimes, such as hillslopes, colluvial channels, and fluvial channels, typically on the basis of variables such as topographic gradient and the upstream contributing area that concentrates runoff (*7*). Whether the flux balance occurs across all these regimes and at all spatial scales (even pointwise) or is only applicable to the total or bulk fluxes at the landscape scale has unavoidable consequences for the dynamic character of the landscape (*8*); the former situation leads to time-invariant (frozen) landforms, whereas the latter allows for a dynamic component of SS landscapes. Although many numerical landscape evolution models result in static SS landscapes under simple boundary conditions (usually vertical uplift and uniform rainfall) (*9*–*14*), physical experiments consistently produce SS landscapes with dynamic landforms (*15*–*18*). This notion of dynamic SS landscapes, where drainage divides continuously migrate and local erosion rates are therefore time-variant and spatially nonuniform, is also supported by field and low-temperature thermochronological evidence (*19*–*21*). Dynamic landscape behavior has been successfully incorporated into some numerical models by various mechanisms, such as landsliding (*22*), the use of more realistic flow-routing algorithms (*23*), or via hillslope-fluvial process interactions (*24*).

If erosion rates vary in space and time, how can one distinguish SS landscapes from transient-state (TS) landscapes, which respond to a change in external forcings? One approach would be to compare the variability in erosion rates of SS landscapes, both in terms of their magnitudes and spatial distribution, with those under TS conditions. Despite good knowledge of how individual landscape components, such as alluvial rivers, bedrock rivers, and hillslopes (*25*–*29*), respond to changes in external forcings, our understanding of the organized erosional response of the landscape as a whole remains elusive. Recent studies have tried to explain the variability of erosion rates in natural landscapes due to, for example, stochasticity of hillslope processes and knickpoint dynamics (*21*, *30*–*32*). However, a comprehensive characterization of this variability, especially in terms of spatial patterns, would demand repeated topographic data at high spatial resolution and over long periods of time. These data are typically not available for natural landscapes, making physical experiments (*15*–*18*, *27*, *33*–*35*) a necessary tool for exploring erosion variability. Although physical experiments have been used to document large-scale TS landscape responses (*15*–*18*, *27*), they have not typically been used to examine the multiscale spatial variability of sediment fluxes under SS conditions to quantify the dynamic nature of SS landscapes and to compare it with TS responses.

Here, we analyze a unique experimental landscape, which provides a detailed space-time record of the topography produced at the eXperimental Landscape Evolution (XLE) facility at St. Anthony Falls Laboratory (see Fig. 1 for the schematic of the XLE facility and Materials and Methods for further description). We seek to (i) fully characterize SS landscapes in terms of local sediment fluxes to advance our understanding of their dynamic nature and (ii) quantify the manner in which landscapes reorganize in response to changes in external forcings.

## RESULTS AND DISCUSSION

### SS landscape

Assuming uniform grain size distribution and material porosity, as is the case in our experiment, the pixel-wise measured topographic change relates to the divergence of sediment flux and the constant uplift rate *U* by the Exner equation(1)

The erosion depth (ED) at pixel *i* over a time interval [*t*, *t* + Δ*t*] is obtained by integrating the flux divergence(2)where positive (negative) values of ED_{i} imply net erosion (deposition) at pixel *i*.

A landscape is said to be at SS when the erosional fluxes balance out the sediment flux provided by the rock uplift. Depending on the scale at which this flux balance is applicable, two different types of SS can be defined (*8*): flux SS and topographic SS. In flux SS, the total flux of sediment leaving the system balances the amount provided by tectonic uplift during an interval of time Δ*t*:(3)where 〈⋅〉 denotes spatial average over all pixels *i* and the first equality acknowledges the time-independent average flux. Flux SS is also referred to as statistical SS, acknowledging that several statistical properties of the landscape, such as slope and upstream contributing area probability distributions, sediment discharge, and river network properties, remain constant (*17*, *18*). In topographic SS, the surface elevation does not change over time because the divergence of sediment flux is the same at every point of the landscape and is exactly equal to the uplift rate(4)

Using the XLE facility, we let the landscape evolve under constant uplift rate *U* = 20 mm/hour and constant precipitation rate *P* = 45 mm/hour for 8 hours. SS conditions were inferred by a time-invariant sediment flux rate equal to the uplift rate (*17*). Figure 2 illustrates the SS nature of the landscape by showing the time invariance of two important statistical properties: the slope-area curve (Fig. 2A) and the probability distribution of pixel-wise ED, which also confirms a constant mean ED (Fig. 2B). The slope-area curves were obtained from four consecutive topographies at SS (measured 5 min apart) using the steepest downslope direction to estimate local slope and the D-infinity algorithm (*36*, *37*) to compute upstream contributing areas. Slope-area curves are a useful tool for revealing the scales of geomorphic organization (*7*, *38*–*45*). From changes in the trends of these curves, we can differentiate three process regimes: hillslopes, draining upstream contributing areas that range from 1 to approximately 10 pixels, or up to 2.5 mm^{2}; a colluvial regime corresponding to intermediate upstream contributing areas of 2.5 to 250 mm^{2}; and a fluvial regime corresponding to upstream contributing areas larger than 250 mm^{2}. The specific values are obtained via analysis of slope increments and detection of change of trends, as discussed in the study by Singh *et al*. (*17*). The overlap of consecutive slope-area curves derived from different topographies at SS shows that there was no significant change in these regimes and thus no structural reorganization of the landscape. Note that the higher variability observed for large upstream contributing areas is due to the smaller sample size available for computing the corresponding slope. We also computed probability density functions (PDFs) of the pixel-wise ED, with positive values indicating erosion and negative values indicating deposition, computed by taking the differences of elevation of consecutive topographies measured 5 min apart. From the overlapping distributions and the results of a Kruskal-Wallis test (see Materials and Methods), we conclude that the PDFs are statistically indistinguishable, revealing the statistical SS nature of erosional and depositional processes. We also note that the shape of the PDFs reveals that the landscape is not frozen (that is, it is not a topographic SS); if it were, the PDF would be just a Dirac delta function (single value) centered at the value of the uplift depth (*U* ⋅ Δ*t*, that is, the depth of material provided by the uplift in Δ*t* = 5 min). The observed complex distribution of local ED raises the question about the spatial distribution of the variability in the erosion magnitude. In the next section, we unveil a stationary scale-dependent pattern of erosion for SS landscapes via a spatial analysis of the sediment fluxes.

### Scale-dependent erosional patterns: The E50-area curve

We ask whether there exists a characteristic erosional signature of SS landscapes reflective of their geomorphologic organization. For this, we interrogate the landscape in terms of the pixel-wise erosion (deposition) depth as a function of the pixel location parameterized by the upstream contributing area. Specifically, we compute the PDF of ED for sets of pixels grouped in 100 equal probability area bins according to their upstream drainage area *A*_{i}. We summarize the results of this analysis in a so-called E50-area curve (Fig. 3A), where we estimate the probability that the pixel-wise ED within each drainage area bin exceeds the median ED of the whole landscape. We highlight two main points revealed by the E50-area curves. First, the stationary shape of the curve for fluxes computed at different SS intervals reveals a statistical pattern that is persistent over time; that is, the E50-area curve is a statistical signature of the SS landscape. Second, the curves have a characteristic nonlinear shape that deviates from the trivial horizontal curve (equal to 0.5 for all values of the upstream contributing area), which would be expected under topographic SS. Specifically, the E50-area curve reveals that the regimes of the landscapes characterized by both small (hillslopes) and large (fluvial) contributing areas erode significantly more than the median of the landscape.

It can seem paradoxical to argue that SS landscapes have a time-invariant erosional signature that is nonuniform across different scales, where, for instance, hillslopes are consistently more likely to erode than the rest of the landscape. This erosional pattern also apparently contradicts the possibility of maintaining the statistical properties of an SS landscape, such as invariant total relief and stationary slope-area curves. The missing factor needed to reconcile these ostensible discrepancies is the dynamic character of the landforms at SS. Asserting that hillslopes are more likely to erode is not equivalent to saying that fixed locations in the landscape are more likely to erode because individual pixels can evolve and belong to different geomorphic regimes at different times. A higher erosion rate in the hillslope pixels reduces their elevation over time and hence changes the upstream contributing areas, eventually shifting them into a regime with a lower erosion rate. To illustrate this dynamic nature of the SS topography, we show in Fig. 3B that 40% of the hillslope pixels (that is, pixels with upstream contributing areas of less than 0.5 mm^{2}) drain larger areas after 5 min of landscape evolution under SS conditions (see fig. S1 for alternative values of initial upstream area). This dynamic behavior ensures that erosion rates estimated using sediment fluxes measured at a fixed location over sufficiently long periods will converge to the erosion rate of the whole landscape, as that fixed location visits different regimes of the E50-area curve.

We emphasize that patterns in erosional fluxes, as shown by the E50-area curve, are easily disguised by examining the landscape in a different manner, for example, by random sampling. Figure 3C shows the probability of erosion for pixels contained in random samples of the same size as those used to build the E50-area curve. The stationarity of the probabilities over time for fixed locations is additional evidence supporting the SS nature of the landscape and, by itself, might lead one to conclude that no persistent spatial patterns of erosion are expected once SS is reached. Figure S2 shows the estimation of erosional rates when different sample sizes are considered, depicting a robust behavior of those estimators for sample sizes even smaller than the one used in Fig. 3C.

The existence of time-invariant spatially explicit patterns of erosion in SS landscapes opens questions of how to detect and characterize the response of the landscape to changing external forcing. In the next section, we show that a similar analysis reveals a significantly distinct hierarchical response of a landscape under increased rainfall intensity.

### TS landscape

A TS landscape can be defined as a landscape with nonzero net material flux at the landscape scale. A TS is normally a consequence of abrupt changes in the external forcings that drive landscape evolution, such as rock uplift rate and precipitation. Using our experimental facility, we investigate the landscape reorganization at the onset of the TS that is produced by a fivefold increase in rainfall intensity. Under TS conditions, the amount of sediment leaving the system significantly exceeds the sediment production provided by tectonic uplift:(5)

Note that depends on both *t* and Δ*t*; the disequilibrium expressed in Eq. 5 gradually decays with time (*17*) as the landscape approaches a new SS.

We are interested in comparing the distinct dynamic response of the reorganizing landscape during the onset of TS conditions with the inherent spatial variability in erosion rates within the SS landscape. However, for a meaningful comparison of the sediment fluxes, the two landscapes must first be rendered comparable in terms of the total volume of sediment that is removed. For this, we integrate the SS and TS landscapes over different time intervals, that is, over a longer time interval (*k*Δ*t*) at SS to match the eroded sediment volume produced over an interval Δ*t* under increased precipitation at TS(6)

Acknowledging the SS condition of Eq. 3, the time-rescaling factor *k*, which depends on both *t* and Δ*t*, can be estimated by the volume rescaling factor, that is, as . Focusing our analysis on the first 5 min (that is, Δ*t* = 5 min) after the transition to increased precipitation rate, we found that *k* = 2.6, meaning that an integration time of 13 min (2.6 × 5 min) is needed at SS to dislodge the same total volume of sediment as that on the first 5 min under TS. This ratio decreases as the integration time increases and eventually approaches *k* = 1 at a new SS (because the uplift rate remains the same). During the experimental run, landscape topography was acquired every 5 min, and so, we can only scale the SS landscape by integer values of *k*. By comparing the PDFs of ED corresponding to different values of *k* (see fig. S3), we select *k* = 2 (that is, topographies measured 10 min apart) as the best estimate within the available temporal discretization for the rest of the study.

The spatial patterns of erosion at TS are substantially different from those at SS (Fig. 4). To quantify the distinct distributed response occurring during the onset of the TS, we show the E50-area curves for SS (Δ*t* = 10 min) and for the onset of the TS (Δ*t* = 5 min), as well as the slope-area curve corresponding to the SS, in Fig. 5A. The E50-area curve at TS shows a significant deviation from that at SS within three distinct regions of erosional regime change under increased precipitation: (i) For areas *A*_{i} < 0.75 mm^{2}, there is a large percentage of high-erosion pixels for both SS and TS, but erosion is enhanced during TS compared to SS; (ii) for areas 0.75 mm^{2} < *A*_{i} < 50 mm^{2}, the percentage of high-erosion pixels decreases with upstream drainage area in both SS and TS, but the rate of decrease is larger in TS than SS; (iii) for areas *A*_{i} > 50 mm^{2}, there is a regime shift from downstream-increasing erosion to downstream-decreasing erosion: Erosion increases sharply with *A* for SS, but for TS, the fraction of highly eroding pixels decreases with *A*. Putting these results in the geomorphic context provided by the slope-area curve, we can conclude that, during landscape reorganization in TS, hillslopes undergo accelerated erosion, colluvial and slightly convergent regions experience reduced erosion, and fluvial channels experience a reduction of their channel incision rate (erosion) due to the increase in sediment flux delivered from upstream. These results are compatible with numerical simulations by Tucker and Slingerland (*10*). Also, note that the emergent scales that demarcate these erosional regime transitions match fairly well with the scales of geomorphic process regime transitions from hillslope to colluvial to fluvial obtained from the slope-area curve (*7*, *44*), as illustrated in Fig. 5A. To the best of our knowledge, this is the first time that these erosional regime transitions (revealed by the E50-area curves) and geomorphic process regime transitions (revealed by the slope-area curves) have been explored simultaneously at the landscape scale to detect and interpret reorganization.

This reorganization can be visualized by explicitly positioning on the landscape all pixels that transition from high to low erosion and vice versa during reorganization, relative to the landscape median erosion rate. Figure 5 (B to D) depicts a single drainage basin and shows the parts of the landscape that have changed their erosional behavior during the onset of TS. It is seen that hillslope pixels are the first to respond to the increased precipitation rate, shifting from low to high erosion values (Fig. 5C). In contrast, fluvial channels shift from high to low erosion values, so that incision rates are reduced because of accelerated upstream erosion and sediment supply (Fig. 5D). Although there is no distinction between sediment and bedrock in our experiment, these results resonate with recent models that suggest that sediment fluxes can exert a significant control on river incision rates (*46*–*49*). The top-down reorganization of the landscape, with information flowing from hillslopes to channels, is distinct to the commonly held view of landscape reorganization in response to base-level changes, in which channels lead and hillslopes follow (*47*, *50*–*53*).

## CONCLUSIONS

The question of whether an SS landscape achieves a frozen topography that exhibits no variability in local erosion rates at any scale or achieves a statistical equilibrium within which erosion dynamically and preferentially changes locally while maintaining the large-scale balance of fluxes remains open. Here, we analyzed a densely monitored experimental landscape to present evidence that SS is characterized by a hierarchical pattern of erosion summarized in a new curve called the E50-area curve. This curve quantifies the probability of a location eroding above or below the landscape median as a function of the location’s upstream contributing area. We explained this curve in terms of the internal dynamics of the SS landscape by showing that locations of the landscape switch geomorphic regimes through time (for example, hillslopes erode more than the landscape median, lowering their relative elevation and increasing their upstream contributing area, thus shifting to a new geomorphic regime). We proposed that the E50-area curve is a characteristic signature of SS landscapes that should be reproduced in numerical models. Finally, we showed how the shape of the E50-area curve changes when the landscape is under TS conditions in response to a change in external forcing. How the shape of the E50-area curve evolves as the landscape approaches a new equilibrium in response to its forcing and whether this new equilibrium differs from the original one are open questions currently under experimental and analytical investigation. Extended experimental data will also allow investigation of the variability of the E50-area curve under different external forcings as an emergent property of landscape organization, informing numerical landscape evolution models and providing important information for quantifying the uncertainty of sampled erosional rates in field landscapes.

## MATERIALS AND METHODS

### Description of the experimental setup

The XLE facility (see Fig. 1 for schematic) consists of an erosion box (0.5 × 0.5 × 0.3 m^{3}) with two main controlling variables: (i) uplift rate, adjusted by lowering two opposing sides mimicking mountain uplift, and (ii) rainfall intensity, simulated using 20 ultrafine misting nozzles (droplet size, <10 μm) to achieve approximate spatial uniformity over the box. The rainfall droplet size was small enough to avoid splash disturbances by the drop impact on the landscape surface. The sediment used in the experiment was a homogeneous mixture of fine silica (*D*_{50} = 25 μm), with ~35% water content by volume. The facility was equipped with a high-resolution laser scanner that could obtain the topographic elevation *h*(*x*,*y*,*t*) of the whole surface in 5 s at a spatial resolution of 0.5 mm and a vertical accuracy of better than 0.5 mm. For this experiment, topographic data were acquired every 5 min. We refer to Singh *et al*. (*17*) for a comprehensive discussion of the experimental setup and collected data.

### Statistical analysis

We used the Kruskal-Wallis test (*54*) to compare the empirical PDFs of the pixel-wise ED obtained by differencing consecutive topographies measured 5 min apart. The Kruskal-Wallis test is a rank-based nonparametric test that does not assume a given distribution for the residuals and can be considered a nonparametric counterpart of the one-way analysis of variance (ANOVA). The sample size (number of pixels) of each ED field was *N* = 677,810.

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/3/9/e1701683/DC1

Scale-dependent erosional patterns in SS and TS landscapes

fig. S1. Dynamic landforms at SS.

fig. S2. Estimation of the probability of erosion larger than the landscape median at SS for different sample sizes.

fig. S3. Comparison of the SS and TS landscapes in terms of the aggregate statistics of ED.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is **not** for commercial advantage and provided the original work is properly cited.

## REFERENCES AND NOTES

**Acknowledgments:**We thank J. Mullin, C. Ellis, and L. Reinhardt for helping with the development of the experimental facility (XLE). We also thank B. Dietrich and J.-L. Grimaud for fruitful discussions during the early stages of this work. We also thank the editor K. Hodges, the associate editor P. Bierman, and three anonymous reviewers for their insightful comments, which helped to improve the focus and presentation of our work.

**Funding:**This research was partially supported by the National Center for Earth-Surface Dynamics (NCED) funded by NSF under agreement EAR-0120914 and by NSF grant EAR-1209402 under the Water Sustainability and Climate Program. A.T. acknowledges financial support from the NCED 2 (NSF grant EAR-1246761) postdoctoral fellowship.

**Author contributions:**All authors contributed equally to this work.

**Competing interests:**The authors declare that they have no competing interests.

**Data and materials availability:**All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

- Copyright © 2017 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).